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ABSTRACT 

We have coupled a fast, parametrized star cluster evolution code to a Markov Chain 
Monte Carlo code to determine the distribution of probable initial conditions of ob¬ 
served star clusters, which may serve as a starting point for future A-body calcula¬ 
tions. In this paper we validate our method by applying it to a set of star clusters 
which have been studied in detail numerically with A-body simulations and Monte 
Carlo methods: the Galactic globular clusters M4, 47 Tucanae, NGC 6397, M22, uj 
Centauri, Palomar 14 and Palomar 4, the Galactic open cluster M67, and the M31 
globular cluster Gl. 

For each cluster we derive a distribution of initial conditions that, after evolution up 
to the cluster’s current age, evolves to the currently observed conditions. We find that 
there is a connection between the morphology of the distribution of initial conditions 
and the dynamical age of a cluster and that a degeneracy in the initial half-mass radius 
towards small radii is present for clusters which have undergone a core collapse during 
their evolution. We find that the results of our method are in agreement with A-body 
and Monte Carlo studies for the majority of clusters. We conclude that our method 
is able to find reliable posteriors for the determined initial mass and half-mass radius 
for observed star clusters, and thus forms an suitable starting point for modeling an 
observed cluster’s evolution. 

Key words: galaxies: star clusters: general - methods: numerical. 


1 INTRODUCTION 

What were the initial conditions of the star clusters we ob¬ 
serve today? Answering this question not only requires accu¬ 
rate observations of the current conditions, but also proper 
modeling of star cluster evolution over a large amount of 
time (for a globular cluster typically 12Gyr), taking into 
account a number of physical processes, such as two-body 
relaxation, three- and four body interactions, the stellar evo¬ 
lution of single and binary stars and the effects of a galactic 
tidal field (Ambartsumian 1938; Chandrasekhar 1942; King 


1958 Hut et al. 1992 Henon 1960 Lee & Ostriker 1987). 


The required complexity of the simulation depends on the 
type of cluster one wishes to model and the physics ques¬ 
tion to be answered. In the Milky Way Galaxy we know 157 
globular clusters (Harris 2010), and about 2200 Galactic star 
clusters in the disc. Of the latter, ~2170 are relatively low- 


mass (< 10 4 Mq) systems which we traditionally classify as 
open clusters ( [Dias et ak 2002 version 3.3 - jan/10/2013), 
while ~30 have masses and luminosities comparable to the 
type of clusters that are typically found in studies of exter¬ 
nal galaxies, classified as young massive (> 10 4 Mq) clusters 


(Portegies Zwart, McMillan, & Gieles 2010). The obvious 
differences between these three cluster types lie in the age, 
the number of stars, the (half-mass) density and the amount 
of gas currently present in the cluster. Each factor compli¬ 
cates the simulation of the cluster evolution, if its value is 
high. However, at the basic level, each of these star clusters 
is subject to the same physical processes ( Larsen| 2002 ) and 
will eventually meet the same fate of complete dissolution 
( Aarseth||1973 Baumgardt|2006 ). 
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1.1 Star cluster physics 

Star clusters form in the collapse of giant molecular clouds 
(see e.g. Alves et al.||2001 Molinari et al.||2014 ). After the 
complex phases of cluster formation and early evolution with 
residual gas expulsion, cluster expansion and re-virialisation 


(2014). The A time complexity makes these simulations 


(see e.g. Baumgardt k Kroupa (2007), Pelupessy k Porte- 

gies Zwart ( 

2012), Banerjee k Kroupa (2013) or Longmore 

et al. (2014; 

)), the star cluster is often assumed to evolve as 


a roughly spherically symmetric gas-less system under the 
influence of the following physical processes: 

(i) the dynamics of the stars on both the large scale (two- 
body relaxation) and on the small scale (three- and four- 
body interactions, binary formation and evolution) (Am- 
bartsumian|1938] |Chandrasekhar||1942 |King||1958 ); 


(ii) the stellar evolution of single stars and binaries flint 
et al. 1 1992); 

(iii) the interaction with the tidal field of the galaxy it 
resides in ( Henon|1960 Lee k Ostriker|1987 ). 


Two-body relaxation tends to drive the cluster to the un¬ 
reachable state of equipartition, where the kinetic energy 
of the stars in the cluster is equalized, see e.g. |Khalisi et al.| 
(2007) and the references therein. This has two major effects: 
1) mass segregation: the most massive stars will sink to the 
cluster center, whereas the lower-mass stars populate the 
halo; 2) core-collapse: due to the decrease of kinetic energy 
in the core, gravitational collapse is no longer supported. 
Hence the core will shrink until high enough core densities 
are reached to produce a new source of kinetic energy: bi¬ 


naries (Aarseth 1973). The binaries will halt the process 


of core-collapse and start the core expansion by interacting 
with less massive stars, in which the former become more 
tightly bound (‘hard’) and the latter escape into the halo 
of the cluster (Aarseth 1973). The half-mass radius will in¬ 


crease driven by stellar evolution and hard binaries (Gieles 
et al.| 2010). Since the low-mass components transferred to 


the halo have higher kinetic energy, this process will also 
cause the preferential loss of low-mass stars. The escape of 
stars over the tidal radius leads to a contraction of the clus¬ 
ter, i.e. the decrease of the half-mass radius, at a (roughly) 
constant density ( Henon|[l961 201 1[ Gieles et al.||2011 ). If, 
on the other hand, the cluster has a significant amount of 
primordial binaries (> 10%), either the core-collapse would 
be less deep in the sense that the core radius would not 
decrease as much as it would without the presense of pri¬ 
mordial binaries, or there would be no core-collapse at all, 


see e.g. Baumgardt et al. (2003a). 


1.2 Simulation techniques 

In computational astronomy great progress has been made 
to develop dedicated codes to study the evolution of star 
clusters. These dedicated codes can roughly be divided in 
three groups ( Alexander k Gieles||2012 ). The first group of 
methods are the A-body simulations. Direct A-body sim¬ 
ulations (e.g. |Aarseth||1973| |Makino||1996| iSpurzem 11999) 


have the smallest number of simplifying assumptions. The 
time complexity of these simulations scales with A -2 , where 
A is the number of stars in the simulation, and solutions 
up to arbitrary precision can be obtained, see Portegies 
Zwart k Boekholt (2014) and Boekholt k Portegies Zwart 


computationally expensive. Globular clusters contain tens 
of thousands to several millions of stars at the present day, 
and theoretical motivations point out that they must have 
had an even greater initial number of stars (Aarseth 1973). 
Because of this large initial number of stars a direct IV- 
body simulation of the evolution of a million body globular 
cluster over its entire lifetime is yet to be performed. Direct 
IV-body simulations with fewer stars that are scaled up 
afterwards to match the observed clusters, have been made 
successfully, and also the direct IV-body integration of small 
globular clusters, or of open clusters, has been done with 


success, see e.g. Baumgardt et al.| (2003b), Zonoozi et al 


(2011), Zonoozi et al. (2014) and |Hurley et al. (2005) for 


IV-body simulations of the evolution of Gl, Palomar 14, 
Palomar 4 and M67, respectively. For the less accurate, but 
more efficient, tree-codes the computation time scales with 


IVlog(IV); see Bedorf k Portegies Zwart (2012) for a review 
of the history of the collisional direct and collision-less 
tree-code IV-body methods. 


Other methods, such as Monte Carlo methods (Metropolis 


k Ulam (1949), Freitag k Benz (2001), Giersz (2006) and 


Giersz et al. (2008)) and Fokker Planck models (Kolmogo- 


roff (1931) and Cohn (1979)), scale both with A and make 
some simplifying assumptions such that the computational 
cost decreases, but at the price of the decrease of the accu¬ 
racy. With this second group of methods, the computation 
of the evolution of globular clusters with initially more 
than 10 6 stars has become possible, but can still take a 
substantial amount of time and computer power. Therefore 
the choice of initial conditions is very important and a lot 
of effort is usually put into choosing a plausible set of initial 
conditions that will evolve to a cluster with characteristics 
similar to the observed cluster of interest. See e.g. |Giersz| 
k Heggie (2009) where a number of small scale runs, i.e. 


simulations for clusters with a lower initial number of stars, 
are performed to this end. 

The third group of methods are the gaseous models 
(Larson 1970) and the semi-analytical methods, where the 
computation time is approximately IV-independent. The 
semi-analytical methods make use of parametric equations 
for the evolution of the cluster variables, such as the 
number of stars, A, and the half-mass radius, rhm, which 
are then solved with a numerical integrator. These methods 
are therefore faster again. One such code is the recently 
developed Evolve Me A Cluster of StarS (emacss) code, 
which is based on the flow of energy in a cluster and allows 
one to compute the evolution of a globular cluster over 12 


Gyr in a fraction of a second (Alexander k Gieles 


2012 


Gieles et al. 2014[ Alexander et al. 2014). This code has 


been tested against A-body simulations and has proven 
to be an extremely powerful tool in understanding cluster 
evolution and exploring large regions in initial parameter 
space. 

1.3 Goal of this work 

We aim to determine the most probable sets of initial 
conditions in total cluster mass, M, and half-mass radius, 
rhm, for any observed star cluster and explore possible 
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degeneracies. Finding the initial conditions of individual 


(2003), Baumgardt et al. (2003b), Hurley et al. ([2005j), 

Heggie & Giersz 

(2008), Giersz & Heggie (2009), 

Giersz & 

Heggie (2011), Zonoozi et al. (2011), Heggie & Giersz (2014) 

and Zonoozi et 

al. (2014). These studies used 

elaborate 


simulation techniques to model cluster evolution and found 
suitable initial conditions that, after evolving up to the 
cluster’s age, resembled a number of observables of their 
star cluster of interest. However, due to the versatility of 
their methods, these studies were only able to investigate 
up to several tens of sets of initial conditions iteratively. 
They could not investigate the uniqueness of their sets of 
initial conditions or explore whether there are multiple, 
significantly different, sets of initial conditions which evolve 
to the current observables, i.e. whether there are degenera¬ 
cies. In this study we aim to address this latter point as 
well. We accomplish this by coupling the fast, parametrized 
star cluster evolution code emacss (Alexand er et al.] [2Q14) 
to the Markov Chain Monte Carlo (MCMC) code emcee 
( Foreman-Mackey et al.|[ 2013). 


In this paper we describe and demonstrate our method. We 
validate it by applying it to nine star clusters that have been 
studied to great extent with either 7V-body simulations or 
Monte Carlo methods and by comparing our results to the 
results of those methods. The paper is organized as follows: 
in Section [2] we explain our method and we summarize the 
functionality of the underlying star cluster evolution code. 
In Section [3] we describe our validation strategy and set out 
the relevant parts of the extensive work that has been done 
on the validation clusters by other authors. In Section [4] we 
show our results and we discuss them in detail. Section [5] 
summarizes the paper and discusses future work. 


2 METHOD 


2.1 The parametrized star cluster evolution code 


We evolve the star clusters using the parametrized star 
cluster evolution code Evolve Me A Cluster of StarS 
(emacss) ([Alexander &; Gieles (2012); Gieles et al. (2014); 
[Alexander et al. ( 2014| ). emacss includes a prescription for 
mass segregation, the evolution of the mean stellar mass, 
m — M/N , as the result of stellar evolution, the resulting 
expansion of the cluster and the escape of stars over the 
tidal radius, r t . After a phase of “unbalanced” evolution, in 
which the evolution and escape of stars are the dominant 
drivers of the evolution, a phase of “balanced” evolution 
starts ( Alexander et ah||2014 ). Here it is assumed that the 
core produces the correct amount of energy to sustain the 
two-body relaxation process. Balanced evolution is assumed 
to start after a fixed number of (half-mass) relaxation times. 


A cluster is evolved on a circular orbit, with constant 
orbital velocity, v, at a constant galactocentric radius, Rgc, 
about the galactic center, emacss assumes a logarithmic 
potential, 0, for the galaxy which imposes a static tidal field 


3 

rj = 


the cluster ( |Gieles et al.|2014 Alexander et al.||2014| ): 

= v 2 ln(RQc) (1) 

GMRl c (2) 


2v 2 


in which rj is the Jacobi radius and G is the gravitational 
constant. Note that rj = rt for the type of potential used 
here. 


The code uses a Kroupa initial mass function (Kroupaj2001) 
with a lower mass limit of 0.1 Mq and an optional upper 
mass limit m up . It was tested against iV-body simulations 
for 77i U p — 15 Mq or 100 Mq with an initial mean mass of 


0.64 Mq in both cases (Alexander et al. 2014). In each of 


our simulations we use the upper mass limit of 100 Mq 
and an initial mean mass of 0.64 Mq. emacss furthermore 
offers an indication of core-collapse based on an ‘average’ 
cluster, which is adequate to determine the state of a 
cluster substantially before or after core collapse, although 
is unreliable at times within a factor ~2 of the predicated 
collapse itself. The code does not explicitly include a 
prescription for stellar interactions and binary formation, 
nor does it account for the effects of a primordial binary 
population. This most recent version of emacss also allows 
one to take into account the effects of dynamical friction. 
We have not included this feature in the simulations in 
this paper, because the studies we compare our results to 
have not included this effect either. We will explore the 
effects of dynamical friction in forthcoming work. The de¬ 
tails of emacss are described in references mentioned above. 


The most recent version of emacss (Alexander et al. 
2014) is available on GitHutj^] and in the Astrophysical 
Mult ipurpos e Software Environment (AMUSE, |Port egies 
Zwart et al. fl2009| l)P1 


2.2 emacss-mcmc method 

We define a set of initial conditions of a star cluster as 
the cluster’s initial mass and initial half-mass radius, i.e. 
(Mi, 7*hm,i)- We constrain the initial conditions for an 
observed star cluster from the observed current mass, M 0 b s , 
half-mass radius, rhm,obs, age, r 0 bs, galactocentric radius, 
RGC,obs and orbital velocity, r Q bs- For each observed cluster 
we simulate n clusters with different initial total masses and 
half-mass radii, Mi and rh m ,i respectively, from t = OGyr 
to t = r 0 bs- These clusters all start out with the same initial 
galactocentric radius and initial orbital velocity, equal to 
the currently observed value^] since both these parameters 
do not change throughout the evolution with emacss 
without dynamical friction. After evolving the clusters, we 
compare their final conditions in mass and half-mass radius, 


1 https: / / github.com / emacss / emacss 

2 http://amusecode.org/ 

3 In this section we describe the general use of the two- 
dimensional version of our method. However, since our aim is 
to validate our method, we choose the same values for the age, 
the galactocentric radius and the orbital velocity as the stud¬ 
ies we compare to in our 2D simulations, but in our 5D sim¬ 
ulations we will marginalize over these three parameters in the 
five-dimensional version of our method, see Section|3.1.2| 
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Initialize ensemble of walkers: 

for walker j : x[ ,0 =[ log ( M-), r hm j] 

x/-°[0]et/(log(Mj,log(M ois )+3) 

x/'°[l]et/(0,500) 


Probability function 


Check parameter validity: 

log(M obs )<log(^-)<log(M obs )+3 

0 <^^< 500 

pc 


Evolve the cluster with EMACSS: 

-for t=t obs ; 

- on circular orbit dXR GC =R GCf0bs \ 


- with static tidal field set by R 


GC,obs’ v obs' 


Check if cluster doesn’t 
dissolve at t^obs- 



Propose next step for walker: 

x J i’*=[\og{M*),r*hm ,,•] 


1 


Pick ueu{ 0 , 1 ) : 

if u < Acceptance Function( p{x{’ k ), p (x/)*:) 


j,k +1 /,* 

X ; = X ; 


else: 


j,k +1 j ,k 

X- = X - 


Compare present conditions of 
simulation to observables (w. errors) 
and assign probability to IC. 

p{xl)ccexp[^-(x J f -\i)ir'(x J f -\l)) 


if k = 0 


if MO 



Figure 1 . The schematic representation of our method. One round through this scheme represents one iteration, both the burn-in and 
subsequent chain iteration phase. 


i.e. Mf and 7*h m ,f respectively, to the observed present day 
values and assign a (posterior) probability to each initial 
condition. Note that the zero age of the cluster, i.e. t = 
0 Gyr, is defined in emacss as the time when both the 
residual gas of the giant molecular gas cloud, from which 
the star cluster formed, has escaped the cluster and the 
cluster has reached virial equilibrium. 


For choosing the initial conditions in mass and half¬ 
mass radius, we could have used a grid of n mass and 
half-mass radius pairs, e.g. evenly spaced in both half-mass 
radius and (the logarithm of) mass. However, if one wants 
to determine the initial conditions in more than two dimen¬ 
sions, which we do in our 5D simulations (see Section 3.1.2), 
the grid approach is no longer feasible and one needs to 
sample the initial conditions with a method that efficiently 


probes and properly covers a multidimensional parameter 
space. We therefore use the affine-invariant ensemble 
sampler for Markov Chain Monte Carlo (MCMC) as coded 
up in the emcee code ( Foreman-Mackey et al.||2013| ) based 


on Goodman & Weare (2010). 


MCMC is a procedure to generate a random walk in 
parameter space to obtain an approximation to the 
unknown posterior density distribution function (PDF) 
(Foreman-Mackey et al. 2013 Metropolis et ah] 1953[ 
Hastings 1970). Sampling the PDF starts by initializing the 
walkers accross parameter space according to some prior 
distribution. The walkers then sample parameter space 
according to the specific MCMC algorithm one employs, 
and eventually converge towards those regions of parameter 
space with high posterior probability. Since in most cases 
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one has no a priori knowledge what the PDF looks like, 
the simplest and most un-informative prior - a uniform 
distribution in each parameter/dimension - is used. The 
walkers burn-in to probable regions of parameter space 
that can then be used a s a starting point of the subs equent 
(chain iteration) phase (For eman-Mackey et al.|2013| p] 


Figure [I] presents the method we employ. One round 
through the scheme represents one iteration. For each of the 
observed clusters we sample a two-dimensional parameter 
space, in log(M) and rh m - We use log(M) instead of M, 
because we experienced that the MCMC method is more 
efficient when it covers the large range of several orders 
of magnitude in mass in logarithmic space than in linear 
space. We use an ensemble of n w walkers, a burn-in phase 
of rib iterations and n c subsequent (chain) iterations such 
that evolve a total of n — n w (l +rib + n c ) clusters. A walker 
j at iteration k is defined as xf k = [log(M/’ fe ), J; the 
subscript i denotes that it concerns an initial condition. We 
use the following boundary conditions for the parameters: 

log(Mobs) ^ lo s(]iF < log(Mobs) + 3, 

0 C < 500, (3) 

in which the lower mass boundary comes from the fact that 
the initial cluster must have been at least as massive as it is 
today. The other boundaries were found to be reasonable val¬ 
ues to not exclude any possibly interesting regions in mass 
and half-mass radius and to get a good balance between 
proper coverage and quick convergence. For the initializa¬ 
tion we have experimented with several different prior dis¬ 
tributions, see Section [4.1 1 When the walkers are initialized, 
they are evaluated in the probability function, see Figure [l] 
The probability function determines how suitable the sets of 
initial conditions are for the observed cluster of interest by 
assigning a posterior probability, p, to each initial condition; 
p takes value in the range 0 to 1 where 0 denotes the low¬ 
est probability and 1 denotes the highest probability. This 
is done as follows: 


(i) Initial condition check: for a walker j at iteration 
k it checks whether its mass and half-mass radius, stored in 
x{' k [ 0] and x{' k [l] respectively, is within the ranges of Eq. 
<§• If a certain initial condition is in the requested range, it 
proceeds to step (ii). If this is not the case, steps (ii)-(iv) will 
be skipped and this initial condition is assigned a posterior 
probability p — 0. 

(ii) Evolution: a cluster with this walker’s initial mass 
and half-mass radius will be evolved with emacss, at a 
constant galactocentric radius .Rgc = ^GC,obs and veloc¬ 
ity v = v 0 bs from t = 0 Gyr until the observed cluster’s age 
t = T 0 bs is reached. It proceeds to step (iii). 

(iii) Dissolution check: it checks whether the cluster 


4 The burn-in is defined as the first part of the Markov chain 
where the walkers reach a certain level of posterior probability, 
but is otherwise the same as the subsequent chain iteration phase. 
Since the first part of the Markov chain usually contains lower 
probabilities, it is common practice to remove the data of the 
burn-in phase from the results such that further analysis is not 
biased by the low probabilities of the burn-in. See e.g. jPntze et al.| 
(2009) for a nice explanation of the burn-in phase. 


dissolved before reaching t = r 0 b s . If for a certain set of ini¬ 
tial conditions the cluster dissolved before reaching the age 
t = T 0 bs, step (iv) will be skipped and this initial condition is 
assigned a posterior probability p = 0. If the cluster stayed 
bound until t = r 0 bs, it proceeds to step (iv). 

(iv) Posterior probability calculation: the final con¬ 
ditions of the cluster - for a walker j at iteration k con¬ 
tained in x j f ,k = [log(Mj’ k ),r^ f ], where the subscript / 
denotes that it concerns a final condition - are compared to 
the observables and a posterior probability p is calculated 
according to: 

ln(p) = -^((4’ k -M) Ts_1 (4’ k -M)). ( 4 ) 

in which ft = [M 0 b s , nim,obs] are the present day observed 
values and E _1 is the inverse of the covariance matrix which 
contains the errors in M 0 b s and rh m ,obs- We assume that the 
errors in mass and half-mass radius are not correlated such 
that the covariance matrix contains no nonzero off-diagonal 
elements: 

_ f\ lo§(Af 0 bs) log(_M 0 bs A Mobs) I 0 \ /r\ 

V 0 A7w,obJ’ 

We chose to use 10% errors in both observables for all our 
clusters in the majority of our simulation; this choice is 
arbitrary and used in this paper as a proof of concept. In 
further applications of our method one can use the actual 
observational errors. In Section |4.3.4| we investigate what 
the effect of taking smaller (5%) or larger (20%) errors has 
on the results. 


After an iteration k each walker will be proposed a 
new set of initial conditions for the next iteration k + 1. 
Whether the walker accepts or declines this proposed 
set is determined as follows: for a walker j a new set of 
initial conditions xf* will be proposed according to the 
stretch move algorithm. This is an algorithm in which one 
simultaneously evolves an ensemble of walkers, whereby 
the proposal distribution, from which the proposed initial 
condition is drawn, is based on the initial conditions for 
the other n w — 1 walkers in the previous iteration (/c), see 
|Foreman-Mackey et al.||~2013| This set of initial conditions 
is evaluated in the probability function, as explained above. 
After this set of initial conditions is assigned a posterior 
probability, p *, p* will be compared to the posterior 
probability of the previous iteration, p , and a probability 
q will be calculated for the acceptance of this proposed 


set of initial conditions, see Eq. (9) of Foreman-Mackey 


et al. (2013). The values of q are in the range 0 to 1 and 


will approach the value Iff p* p and the value 0 if 
p* <C p. Lastly, a random number u will be drawn from a 
uniform distribution between 0 and 1; if q > u the walker 
j accepts the proposed set of initial conditions such that 
x{' k+1 — x{’*, i.e. the walker ‘walks’ to this proposed set. 
If q < u, then x{' k+1 — x\' k and the walker ‘stays’ at its 
previous set of initial conditions. This procedure is repeated 
for all the iterations, both the burn-in iterations and the 
subsequent chain iterations. See |Foreman-Mackey et al.| 
(2013) for further details of the code emcee. 


In Section 14.11 we first determine suitable values for 
the number of walkers, n w , the number of burn-in itera¬ 
tions, rib, and the number of subsequent chain iterations, 
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n c , while testing the performance of our method. Thereafter 
we investigate the effect of different prior distributions on 
the determined distribution of probable initial conditions. 


2.3 Correcting observations for mass segregation 

When one wishes to compare observed structural param¬ 
eters to the simulated ones, it is of great importance that 
these parameters are obtained in the same way. Since it is 
observationally more practical to find the size of a cluster 
by determining the angular projected half-light radius in 
arcminutes (arcmin), whereas theoretically half-mass radii 
in parsec (pc) are more functional, some conversions are 
needed before a proper comparison between the simulations 
and the observations can be made. 


The angular projected half-light radius, 6 P hi, in arcmin 
from the observations can be converted to the projected 
half-light radius, r p hi, in pc once the distance from the 
Earth to the cluster, Tfe, is known: 

rphl = (pc ) tan (aranin ' Moo)' (6) 

We calculate the observational half-mass radius, rhm, for 
two extreme cases: with or without a correction for mass 
segregation (MS). 

(i) Case 1: no correction for MS 
We assume that the cluster experienced no MS at t — r 0 bs 
yet such that the 3D half-mass radius, rhm, is equal to the 3D 
half-light radius, r*hi. In this case we convert the projected 
half-light radius to the 3D half-light radius by multiplying it 
with a geometrical factor 4/3, correcting for the projection: 


^hm — ^hl — (4/3)r p hl 


(7) 


(ii) Case 2: with a correction for MS 
We assume that the cluster did experience an amount of 
MS at t — T 0 bs in the form of a constant conversion factor 
between the projected half-light radius, r p hi, and the 3D 
half-mass radius, rhm, which we read off from Figure 4 of 


Hurley (2007): 


rhm — Crns^phl, 


( 8 ) 


with c ms = 1.9 for clusters with ages > 7 Gyr and c ms = 
1.8 for a cluster ~ 4 Gyr of age. Note that this conversion 
includes the geometric conversion factor of 4/3 and a factor 
~ 1.425 respectively ~ 1.35 to account for MS. 


By doing this we can compare the observational half-mass 
radii to the half-mass radii from the simulations. See the 
third and fourth column of Table [2] for the calculated obser¬ 
vational half-mass radii without and with a correction for 
mass segregation, respectively. 


2.4 Confidence regions 

After a simulation we obtain sets of initial conditions, final 
conditions and their corresponding posterior probabilities 
for n clusters. This also includes proposed sets of initial 
conditions, even if these were eventually rejected. From 
these n initial conditions we remove the initial conditions 
from the initialization, from the burn-in and those outside 


the ranges given in Eq. §• Hence we have each partic¬ 
ular initial condition appear only once in our data. The 
remaining initial conditions include those which survive 
until t = T 0 bs and those which dissolve before reaching the 
observed cluster’s age. Besides analyzing the most probable 
regions in initial total mass versus initial half-mass radius, 
we namely also want to study the regions which are not 
suitable for producing the currently observed clusters. 

Of the surviving initial conditions, we determine the 
regions with a 68.3% and 99.7% confidence level: 1) we 
calculate the normalized posterior probability of each of 
the sets of initial conditions by dividing the posterior 
probability of each set by the sum of all the posterior 
probabilities; 2) from the set of initial conditions with 
the highest normalized posterior probability downwards, 
we sum the normalized posterior probabilities of each 
subsequent set of initial conditions until this sum equals 
0.997 (0.683). The sets of initial conditions included in that 
sum are a subset of initial conditions which make up 99.7% 
(68.3%) of the total posterior probability, i.e. the 99.7% 
(68.3%) confidence region. Note that our definition of the 
99.7% confidence region is different from what is usually 
meant with a 99.7% confidence region, namely the region 
containing 99.7% of all the data points. The reason for 
our alternative definition is that our aim is to show those 
regions of initial conditions with high posterior probability 
and if we were to use 99.7% of all the data points, this 
region would contain a large number of initial conditions 
with low posterior probability (p < 0.01). 

In Figure [2] we show a simple example of our results 
for the cluster M4 without a correction for mass segre¬ 
gation, whereby the 99.7% confidence region in both the 
initial and the final conditions are enclosed by the yellow, 
dotted contours. 


3 VALIDATION 

3.1 Validation strategy 

We validate our method by applying it to nine star clusters 
that have been studied to great extent with either iV-body 
simulations or Monte Carlo methods. These nine clusters 
include one Galactic open cluster, seven Galactic globular 
clusters and one extragalactic globular cluster; see Table [l] 
for the names of the clusters, the cluster types, their host 
galaxy, the references to papers in which these clusters were 
studied and with which simulation technique. 


3.1.1 Direct comparison 

The first step in the validation is to do a direct comparison of 
the cluster evolution codes by starting from the same sets of 
initial conditions and evolving the cluster under conditions 
which are as similar as possible. Therefore, for each cluster 
we start at the (best-fit) set of initial condition put forward 
by the dedicated studies and we evolve this cluster up to 
the same age under similar conditions, which we describe in 
Section T3. 2 1 and summarize in Table |2j 
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# lo g( r hm/P C ) iogtam/pc) # 


Figure 2. A simple example of the results for the cluster M4. Initial and final condition distributions in half-mass radius versus total 
mass, where the three panels on the left show the initial conditions and the three panels on the right show the final conditions. The 
square panels show the two-dimensional histograms in log(M) vs log(rh m ) of the 99.7% confidence region of the simulations that fit the 
observables without a correction for mass segregation in a black-white density plot: the darker the area, the more initial conditions were 
sampled in this area. Over-plotted on these two-dimensional histograms are the outer contours of the two-dimensional histograms of the 
99.7% (yellow dotted line) confidence level initial conditions for both the simulations without a correction for mass segregation, with the 
projected histograms for log(rh m ) at the top and for log(M) on the left respectively right side. The yellow square with error bars show 
the observation of the cluster’s current mass and half-mass radius when no correction for mass segregation is made, and the cyan square 
with error bars shows the observable with a correction. The pink filled circle in the left panel denotes the initial condition used by the 
dedicated study, which evolves to the final condition shown by the pink filled circle in the right panel. 


Name 

Cluster type 

Galaxy 

Simulation technique 

Reference 


M67 

open cluster 

Milky Way 

Direct A-body 

Hurley et al. (2005 ) 

NGC 6397 

globular cluster 

Milky Way 

Monte Carlo 

Ciersz & Heggie 

(2009 

) 

M4 

globular cluster 

Milky Way 

Monte Carlo 

Heggie &; Ciersz 

(2008 

) 

M22 

globular cluster 

Milky Way 

Monte Carlo 

Heggie & Ciersz 

(2014 

) 

Palomar 4 

globular cluster 

Milky Way 

Direct A-body 

Zonoozi et al. (2014) 

47 Tuc 

globular cluster 

Milky Way 

Monte Carlo 

Ciersz & Heggie 

(2011 

) 

Palomar 14 

globular cluster 

Milky Way 

Direct A-body 

Zonoozi et al. (2011) 

uj Cen 

globular cluster 

Milky Way 

Monte Carlo 

Ciersz & Heggie 

(2003 

1 

G1 

globular cluster 

M31 

Scaled TV-body 

Baumgardt et al.f(200^ 

ibi 


Table 1 . The nine star clusters on which we apply our method to validate it. The first column lists the name of the cluster, the second 
column the cluster type and the third column the galaxy in which the cluster resides. The fifth column shows the references to the studies 
which already studied these clusters to great extent with the simulation technique mentioned in column four. The clusters are mentioned 
in the order of increasing current half-mass relaxation time that was taken from the Harris Catalogue ( |Harris ]2010), except for the first 
and last mentioned cluster, where we took the relaxation times from the corresponding reference in this table. See Table [2]for the values 
of these relaxation times. 


3.1.2 Independent emacss-mcmc runs 

The second step is to run our emacss-mcmc method (Sec¬ 
tion 2.2) independently and determine which initial masses 


are presently better constrained by more recent studies, see 


and half-mass radii reproduce the cluster observables best, 
explore possible degeneracies therein and observe whether 
the best-fit initial conditions of the dedicated study are 
contained within our confidence regions. For the sake of 
comparison to the mentioned dedicated studies, we adopt 
the cluster observables that were used in the references 
mentioned in Table [l] even though some cluster parameters 


e.g. Table 1 of the recent study of Marks & Kroupa (2010), 
where present day masses and half-mass radii are listed for 
a number of clusters. Our goal in this paper is to validate 
our method by comparing our results to the results which 
were obtained with the dedicated studies. Therefore it is 
more important to adopt the observables from these studies 
than to use the currently best constrained observables. See 
Section 13.21 and Tabled 
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Name 

Mobs 

[10 4 M 0 ] 

rhm,obs 

[pc] 

^hm,obs (MS corr.) 

[pc] 

-RGC.obs 

[kpc] 

T obs 

[Gyr] 

r'o b.s 

[km/s] 

log(irh) 

[yr] 

M67 

0.14 

3.35 

4.52 

8.0 

4 

220 

8.48 

NGC 6397 

6.6 

2.30 

3.28 

1.99 

12 

220 

8.6 

M4 

6.3 

2.43 

3.47 

1.68 

12 

220 

8.93 

M22 

33 

4.17 

5.94 

5.28 

12 

220 

9.23 

Pal 4 

2.98 

24.7 

35.2 

102.8 

11 

200 

9.42 

47 Tuc 

110 

4.87 

6.94 

2.95 

12 

220 

9.55 

Pal 14 

1.1725 

35.4 

50.2 

71.6 

11 

220 

10.02 

uj Cen 

390 

9.89 

14.09 

1.22 

12 

220 

10.09 

Gl 

1500 

17.6 

25.08 

- 

12 

- 

10.7 


Table 2. Observational data of the validation clusters as used in the references mentioned in Table [I] The first column lists the cluster 
name and the second column the observed mass in 10 4 Mq. The third and fourth columns lists the half-mass radius without a correction 
for mass segregation (i.e. equal to the 3D half-light radius) and the half-mass radius with a correction for mass segregation (MS) 
respectively. The galactocentric radius, age, orbital velocity of each cluster as used in the references mentioned in Table [l] are listed in 
column five, six and seven respectively. For the sake of comparison to these dedicated studies, we use these observables in our simulations, 
even though some cluster parameters are presently better constrained by more recent observations. See Section |3.2| for the reasoning 
behind choosing each of these parameter values. 


The next step is to exploit the power of MCMC. In 
the five-dimensional version of our emacss-mcmc method 
we add the galactocentric distance, age and orbital velocity 
as nuisance parameters and determine the most probable 
initial masses and half-mass radii. By comparing this to the 
results of the two-dimensional runs explained above we can 
explore how dependent probable initial masses and half¬ 
mass radii are on the choice of the galactocentric radius, 
Rgc, age, t, and velocity, v. Adding these parameters as 
nuisance parameters is observationally motivated, since all 
observables are always determined with some error. The 
method here is similar to the description in Section |2.2| 
but the difference is that now at each iteration we also 
sample a galactocentric radius, an age and an orbital 
velocity from Gaussian distributions with the mean equal 
to i?GC,obs, T 0 bs and v Q b s respectively, as adopted from the 
references mentioned in Table [l] see Table [2] and a standard 
deviation equal to 10% of the mean. For Gl, though, we 
take i?GC,obs = 40kpc as mentioned in Baumga rdt et aL\ 
(2003b) and v 0 b s = 230km/s. We perform our 5D runs only 
without a correction for mass segregation. 


The final step is to investigate the stability of our deter¬ 
mined initial conditions against errors in the observed data. 
We do this in two ways. We first compare our 2D. The sec¬ 
ond thing we do, is checking how increasing/decreasing the 
errorbars on the fitting parameters (log(M) and rh m ) will 
effect the distribution of initial conditions. The second thing 
we do, is testing the stability of the initial condtions against 
varying the parameters Rgc, v and r within their errorbars, 
assuming that the observed mass and half-mass radius are 
unchanged. We take 10% errors and look at the maximum 
difference, e.g. comparing a simulation with Rgc = Rg c,obs 
to simulations with Rgc = 7tMc,obs + O.lifcc^bs and 
Rgc — #GC,obs - 0. l-R G c,obs respectively. 


3.2 Validation clusters 

3.2.1 uj Centauri, M4, NGC 6397, 47 Tucanae and M22 


In this section we describe the studies using the Monte 


Carlo code MOCCA ( 

Giersz 1998 2001 

2006; 

) modeling five 

Galactic globular clusters: uj Cen (NGC 5139; 

|Giersz & Heg- 


gie|2003), M4 (NGC 61 21; |Heggie fe Giersz|2008), NGC 6397 
( Giersz fc Heggie|2009| l, 47 Tuc (NGC 104; Giersz & Heggie 


2011) and M22 (NGC 6656; |Heggie fe Giersz||2014) . Her& 


after we refer to these five studies as GH03-14. 


3.2.1.1 Simulation Technique For each of the globu¬ 
lar clusters GH03-14 determined a set of initial conditions 
- and a few sets of initial conditions in the case of M22 - 
by means of small scale modeling, i.e. by performing simula¬ 
tions for clusters with a lower initial number of stars. These 
initial conditions were subsequently evolved up to the age of 
the cluster, which was taken to be 12 Gyr for each of these 
clusters. For M4, NGC 6397, 47 Tuc and M22 their simu¬ 
lations included prescriptions for single and binary stellar 
evolution, for the Galactic tidal field, for two-body relax¬ 
ation and for binaries and their dynamical interactions. For 
uj Cen, which was studied with an early version of MOCCA, 
their simulation included simple prescriptions for single stel¬ 
lar evolution and for the dynamical evolution in the Galactic 
tidal field. After evolution these clusters gave a satisfactory 
match to a number of observed characteristics, such as the 
surface brightness profile and the velocity dispersion profile. 
See the references mentioned in Section [3.2.11 for the details 
of their MOCCA code. 


3.2.1.2 Orbit and tides GH03-14 evolved each of 
these clusters on a circular orbit with a circular velocity of 
220 km/^] They set the tides by imposing an initial tidal 
radius, r*t,i. Given this tidal radius and their initial total 
cluster mass, M, and assuming an isothermal model for the 
Galaxy, we calculate per cluster with which galactocentric 
radius, Rgc, their model is consistent by combining Eq.s 
0 and 0 - We use these galactocentric radii as input for 
our simulations, see the fifth column of Table [2j 

They evolved their model clusters with an initial tidal 
radius of 35 pc, 86 pc, 40 pc, 89 pc and 90 pc, which are 
consistent with a galactocentric radius of 1.68 kpc, 2.95 kpc, 
1.99 kpc, 5.28 kpc and 1.22 kpc for the clusters M4, 47 Tuc, 


5 |Giersz &; Heg gie (2003) do not mention which circular velocity 
they adopted, but we assume it was a value of 220 km/s as well. 
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NGC 6397, M22 (we compare to their model B, see Table 1 
of Heggie & Giersz ( 2014| )) and uj Cen, respectively. Given 
the fact that M4 is on an eccentric orbit with (eccentricity 
e, perigalacticon R p / kpc, apogalacticon R a /kpc) of about 
(0.8, 0.6, 5.9) ( |Dinescu et al. |l999), this cluster experiences 
strong tides near perigalacticon. Therefore their choice of 
the initial tidal radius, and its corresponding galactocentric 
radius, are very reasonable and better than evolving the 


cluster at its current galactocentric radius of 5.9kpc (Harris 


2010): the averaged tidal field experienced by their cluster 


during its 12 Gyr circular orbit with a galactocentric radius 
of 1.68 kpc is comparable to the tidal field experienced by 
M4 in its actual eccentric orbit. Their choices for the initial 
tidal radii for NGC 6397 and M22 with (e, R p , R a ) of about 
(0.34, 3.1, 6.3) and (0.53, 2.9, 9.3) respectively ( |Dinescu 


|et al.[ p~999) are also reasonable, but their choices for uj 
Cen and 47 Tuc with (e, R p , R a ) of about (0.67, 1.2, 6.2) 
and (0.17, 5.2, 7.3) respectively, are slightly overestimating 
the effect of the tidal field. However, in order to have a 
good comparison to the studies mentioned above, for each 
of these clusters we evolve all our model clusters on a 
circular orbit with a circular velocity of 220 km/s at the 
galactocentric radius to which these models are consistent 
with, see Table [2] 

3.2.1.3 Observed mass and radius GH 03-14 get t he 
observational mass of M4 of 6.3 • 10 4 Mq from 


(2004) and got the observed half-light radius from 


Richer et al. 


Harris 


(1996), which in turn lists angular half-mass radii, stated to 


be taken from the direct average of Trager et al. (1995) and 


van den Bergh et al. (1991). Since both latter studies obtain 


their half-light radii from projected data, we assume that 

on between half-light radii 
(2010) lists half-light radii 


Harris (1996) used a 1-to-l relation between half-light radii 
and half-mass radii, since 


Harris 


and these values are in many cases very similar to the in 


Harris (1996) mentioned ‘half-mass radii’. The 2.3 pc for the 


half-light radius mentioned in Table 2 of Heggie ^fe Giersz| 
(2008) is a mistyped value (Heggie, private communication), 
since the angular projected half-light radius 0 P hi given in 


Harris 


(1996) is 3.65 arcmin (see also Table 1 of Giersz & 


Heggie (2009)). In combination with the distance between 
M4 and the Earth, Re = 1.72 kpc mentioned in |Heggie 
Giersz (2008), 0 p h\ converts to r p hi = 1.83 pc using Eq. |6j). 
We therefore use r p hi = 1.83 pc; for the conversion to the 
3D half-mass radius rh m with a correction for mass segrega¬ 


tion (see Section 2.3), we use the relation (18} with cms = 1.9. 


GH03-14 got th e observation al mass for 47 Tuc of 
1.1 • 10 6 Mq from |Meylan| < |l989| ) , for NGC 6397 of 6.6 • 10 4 
M 0 from Drukierj! 19951) and for uj Cen of 3.9 • 10® M 0 


from Pryor & Meylan (|1993). They do not mention the 


observational cluster mass of M22, so we take the cluster 


mass of 3.3 • 10 5 M 0 from Richer et al. (2008). GH03-14 


got the observed radius for 47 Tuc of 2.79 arcmin and for 


NGC 6397 of 2.33 arcmin both from Harris (1996); they do 


not mention the observational radius for uj Cen and M22, 
so we get the observed radius for uj Cen of 5.00 arcmin 
from 


Harris 


(1996) and for M22 of 3.36 arcmin from 


Harris 


(2010) respectively. Using the same reasoning as mentioned 
for M4, we assume the mentioned observed radii to be the 
projected (2D) angular half-light radii 0 p hi. We calculate 
r ph i in by using Eq. § in combination with distances of 


47 Tuc, NGC 6397, uj Cen and M22 to the Earth, Re, of 
4.5kpc ( Giersz &; Heggie||2011 ), 2.55kpc (Giersz & Heggie 
2009), 5.1 kpc (|Harris |1996} and 3.2 kpc (|Harris| |2010} 


respectively. For the conversion to the 3D half-mass radius 
rhm with a correction for mass segregation, we use the 
relation © with cms — 1.9. 

3.2.1.4 Model initial and final conditions The best- 
fit sets of initial conditions GH03-14 found for NGC 6397, 
M4, 47 Tuc and M22 (model B) in (mass/M 0 , half-mass 
radius/pc) are listed in Table [3] For uj Cen Giersz & Heg¬ 


gie] (2003) mention the initial and final mass of their best-fit 
model (we list those in Table[3]), but they do not mention the 
initial and final half-mass radii of that model, only their ini¬ 
tial and final tidal radii. Assuming that the tidal radius is the 
edge radius of the King model initially (Heggie, private com¬ 
munication) , we can calculate the initial half-mass radius by 
using the ratio of tidal-to-half-mass radius taken from Fig¬ 
ure 8.3. in |Heggie~fc Hut (2003). Using r t /rh m ~ 9.65 for a 
central potential Wo = 7.7, we find rh m ,i = 9.33 pc. It is not 
clear whether this assumption is still valid after 12 Gyr of 
evolution, so we cannot calculate their final half-mass radius. 


3.2.2 Palomar 14 and Palomar 4 

In this section we describe the first published direct iV-body 
simulations of the two large and sparse Galactic globular 
clusters, both residing in the outer halo: Pal 14 (Zonoozi 
et al.pUTT ) and Pal 4 ( Zonoozi et al.||2014 ). They use the 


collisional iV-body code NBODY6 (Aarseth 


2003) on GPU 


computers. Hereafter we refer to these two studies by Zll- 
14. 

3.2.2.1 Simulation Technique Zll-14 simulate the 
evolution of Palomar 14 and Palomar 4 and for Pal 14 they 
compute 65 models and for Pal 4 a total of 20 models, di¬ 
vided in three categories: 1) clusters with a Kroupa (2001) 
IMF in the range 0.08 < M/M 0 < 100 (referred to as their 
canonical-NS model), 2) clusters with a flattened IMF, 3) 
clusters with a [Kroupa (2001) IMF, but with primordial 
mass segregation. For Pal 14 they computed one additional 
model with a Kroupa (2001) IMF, but with primordial bi¬ 


naries. See the references mentioned in Section 13.2.21 for the 
details of their simulations. 

3.2.2.2 Orbit and tides For Pal 14 Zll-14 varied 
the initial half-mass radius and mass and evolved each 
cluster for 11 Gyr on a circular orbit in a logarithmic 
potential with a circular velocity of 220 km/s at Palomar 
14’s currently observed galactocentric radius, of which the 
authors do not mention the value they adopted. After the 
evolution, they fit the final values for the number of bright 
stars, the projected half-light radius and the slope of the 
mass function in the mass range 0.525-0.795 M 0 to the 
observed ones. To make a fair comparison, we also evolve 
all of our model clusters for Pal 14 on a circular orbit at its 
current galactocentric radius, for which we adopt a value 


of 71.6kpc from Harris (2010) with a circular velocity of 


220 km/s. To the best of our knowledge there are no refer¬ 
ences reporting on the orbit of Pal 14. Jordi et al. (2009) 


mention that the orbit could possibly be eccentric, in which 























































































































































10 J. T. Pijloo et al. 


case evolving the cluster on a circular orbit at the current 
galactocentric radius is an underestimation of the tidal field. 

For Pal 4 Zll-14 also varied the initial half-mass ra¬ 
dius and mass and evolve each cluster for 11 Gyr on 
a circular orbit with a circular velocity of 200 km/s at 
Palomar 4’s current galactocentric radius of 102.8 kpc in 
an analytic galactic background potential consisting of a 
bulge, a disc and a logarithmic halo, which they adjusted 
to resemble the Milky Way. 

3.2.2.3 Observed mass and radius Zll-14 got the 

observed total mass for Pal 14 of about 12000 M© fro m|Jordi 


half-mass relaxation time t r h for their model as was inferred 


et al. 


et al. 


2009) and for Pal 4 of about 29800 M© from Frank 


2012). For Pal 14 Zll-14 got the observed projected 


angular half-light radius 0 P hm of 1.28 arcmin from |Hilker| 
(2006), which they in turn convert to the projected (2D) 


half-light radius r p hi of 26.4 ± 0.5 pc and a 3D half-light 
radius rhi of 35.4 =b 0.6 pc. For Pal 4 Zll-14 got the observed 
projected angular half-light radius 6 P h m of 0.62 arcmin from 
King (1966) model fitting on WFPC2 data and broad-band 


imaging with the Low-Resolution Imaging Spectrometer at 
the Keck II telescope, which they convert to the projected 
(2D) half-light radius r p hi of 18.4 =b 1.1 pc and a 3D half-light 
radius rhi of about 24 pc. We converted these projected half- 
light radii to the 3D half-mass radii rh m with a correction 


for mass segregation (see Section 2.3) by using the relation 
© with cms — 1.9. 

3.2.2.4 Model initial and final conditions For both 
Pal 14 and Pal 4 we compare to their canonical-NS models, 
since these models are most comparable to emacss. For Pal 
14 these are the 28 models mentioned in Table 1 of lZonoozH 
|et al. ([201l]) with initial masses in the range 40000-60000 M© 
(Zonoozi, private communication) and initial half-mass radii 
in the range 15-25 pc. For Pal 4 these are the seven models 
mentioned in Table 1 of |Zonoozi et ah (2014) with initial 
masses in the range 50000-57000 M© and initial half-mass 
radii in the range 12-14.5 pc (Zonoozi, private communi¬ 
cation). The fact that these are not their best-fit models 
is not a problem for validation purposes. However, in our 
Figure [lO] and [T2j showing the results of the independent 
emacss-mcmc runs in 2D, we also plot all the other models 
of Zonoozi et al.| ( |2011| and |Zonoozi et aL (2014) to see if 
their other (including best fitting) models are contained in 
our confidence regions. 


3.2.3 G1 

In this section we describe the study using scaled iV-body 
modeling to investigate the evolution of M31’s largest glob¬ 
ular clusters G1 (Mayall II; Baumgardt et al.|2003b ). They 
use the collisional iV-body code NBODY4 (Aarseth 1999) 
on the GRAPE-6 computers ( Makino et al.||2003 ). 


3.2.3.1 Simulation Technique |Baumgardt et al.|2003b| 
simulate the evolution of G1 by running simulations for star 
clusters with N = 65536 stars, since direct iV-body sim¬ 
ulations with a number of stars similar to th e number of 
stars present in G1 (N ~ 10 7 according to Baumgardt et al. 


for G1 from observations (Meylan et al. (2001) estimated 
t r h ~ 50 Gyr). They perform several dozen of runs to de¬ 
termine their best fit to the surface density, velocity dis¬ 
persion, rotation and ellipticity profiles. Baumga rdt et ah'] 
2003b| constructed two models: 1) a single non-rotating clus¬ 
ter, and 2) a rotating merger product, where the merger 
occurred during the formation process. They varied the ini¬ 
tial density profiles, half-mass radii, total masses, and global 
mass-to-light ratios M/L and evolve each cluster for 13 Gyr. 
The authors construct the final density and velocity pro¬ 
files from 10 snapshots in the range 11.75 - 12.25 Gyr and 
give their final cluster mass and half-mass radius at 12 


Gyr. They use a Kroupa (2001) mass function in the range 


0.1 < M/M© < 30 and include the effects of stellar evo¬ 
lution and two-body relaxation. Their simulations do not 
contain primordial binaries, as G1 should still be far from 
core-collapse. See the reference mentioned in Section |3.2.3| 
for the details of their simulation. 


3.2.3.2 Orbit and tides [Baumgardt et al.|2003b| do not 
include a tidal field, since they argue that the tides would 
have a negligible effect on the cluster’s evolution, since the 
cluster is currently at a distance of 40 kpc to the center of 
M31 (Meylan et al. 200T|) . We therefore evolve all our model 
clusters for G1 for 12 Gyr without including a tidal field (i.e. 
as an isolated cluster); see Table [ 2 ] 


3.2.3.3 Observed mass and radius Baumgardt et al. 
(2003b) got a set of observed half-mass radii rhm in the 
range 12.3 - 15.0 pc and observed total masse s in the range 

7.3 - 17 -10 6 M© from Meylan et al. (2001), who in turn 


estimated the half-mass radii and masses from the surface 
brightness profile from HST/WFPC2 images and velocity 
dispersion profile from KECK/HIRES spectra in combina¬ 
tion with King model, King-Michie model and virial theorem 
estimates. We choose to compare our results to the observ¬ 
ables which they obtain with their King-Michie model nr 
4, which gives somewhat average values of the above men¬ 
tioned ranges: rhm = 13.2 pc, M = 15 • 10 6 M© and t r h ~ 50 
Gyr. However, we have checked that the ‘half-mass radius’ 


mentioned in Meylan et al. (2001) comes from their angular 
projected radii in arcmin by using Eq. © in combination 
with its distance to the Earth Re = 770 kpc, which they 
provided in their Table 1. They do not mention that they 
corrected for projection (factor 4/3) or that they did any 
correction for mass segregation. We therefore assume that 
the radius they refer to as the half-mass radius is actually 
the projected half-light radius. For the conversion to the 3D 
half-mass radius rhm with a correction for mass segregation 
(see Section 2.3), we use the relation © with cms — 1-9. 


(2003b)) was, and still is, out of reach. They used the same 


3.2.3.4 Model initial and final conditions For G1 

we compare to their non-rotating model, since this model is 
most comparable to emacss. After scaling up, the cluster of 
this model obtained a final mass of 7.6 TO 6 M© and a final 
half-mass radius of 6.76pc (Baumg ardt et al7||2003b ). They 
found that during the evolution the cluster mainly expanded 
by a factor of 1.75 due to stellar evolution and that it lost 
about 40% of its mass over 12 Gyr (Baumgardt, private 
communication). From this we calculated an initial mass of 
1.27 TO 7 M© and an initial half-mass radius of 3.86 pc. 
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3.2.4 M67 

In this section we describe the work using direct iV-body 
modeling to study the evolution of the rich and relatively 


old Galactic open cluster M67 (NGC 2682; Hurley et ah 


2005 ). The authors simulate the evolution of M67 by using 
the collisional iV-body code NBODY 4 (|Aarseth|1999 ) on the 
GRAPE-6 computers ( Makino et al.||2003 ). 


et al. (2005) estimated that this luminous mass represented 


a total cluster mass of about 1400 M 0 . 

3.2.4.4 Model initial and final conditions The best- 
fit initial conditions Hurley et al. (2005) found for M67 
(model 2) in (mass/M©, half-mass radius/pc) are listed in 
Table H 


3.2.4.1 Simulation Technique Hurley et al. (2005) 


modeled the evolution of M67 by performing iV-body simu¬ 
lations. They compared their modeled surface density pro¬ 
file to the surface density profile of M67 of Bonatto & Bica 
(2005) provided by Bonatto (private communication), their 
modeled color-magnitude diagram (CMD) to the observed 
CMD of Montgomery et al. ( |1993| and their modeled lumi¬ 
nosity function and their structural parameters such as the 
half-mass radius to the observational data from IFan et ahl 


(1996). They furthermore extensively study the stellar pop¬ 


ulations in their simulation and especially focus on the for¬ 
mation channels of blue stragglers (BSs) and compare their 
results to observat ional data from|Fan et al. (1996), Latham 


& Milone (1996), Milone &; Latham| (|1992 ) and Leonard 


(1996). For the single stars, they use a Kroupa et al. (1993) 


mass function in the range 0.1 < M/M© < 50 and their 
model fully accounts for the effects of cluster dynamics as 
well as stellar and binary evolution, including a significant 


fraction of primordial binaries. Hurley et al. (2005) con¬ 


structed two different models, differing only in the initial 
mass, the Galactocentric radius and the binary period dis¬ 
tribution. Their second and favoredcolor model are ran for 


a star cluster with N = 36000 stars. See Hurley et al. (2005) 


for the details of their models and the iV-body code they 
used. 


3.2.4.2 Orbit and tides We compare to the second 


model of Hurley et al. (2005), because it is their best-fit 


model. In this model they evolved the cluster for 4 Gyr on a 
circular orbit with a circular velocity of 220 km/s at a galac¬ 
tocentric radius of 8.0 kpc, which is a reasonable choice for 
a cluster on an slightly eccentric orbit with a perigalaciticon 


of 6.8 kpc and a apogalacticon of 9.1 kpc (Carraro & Chiosi 
1994). 


3.2.4.3 Observed mass and radius The half-mass ra¬ 
dius of Main Sequence stars observed within 10 pc th at|Hur-| 


ley et al. (2005) used was taken from Fan et al. (1996), 


who determined it to be 2.5 pc. However, we checked that 
the ‘half-mass radius’ mentioned in Fan et al. (1996) comes 
from converting their angular projected radii in arcmin by 
using Eq. © in combination with the cluster’s distance to 
the Earth, Re = 783 pc, calculated from their provided dis¬ 
tance modulus of 9.47 mag. They do not mention that they 
corrected for projection (factor 4/3) or that they did any 
correction for mass segregation. We therefore assume that 
the radius they refer to as the half-mass radius, is actually 
the projected half-light radius. For the conversion to the 3D 
half-mass radius rh m with a correction for mass segregation 
(see Section 2.3), we use the relation © with cms — 1.8. 
Both studies took the total luminous cluster mass from lFanll 


et al. (1996), which determined that to be 1000 M©. Hurley 


4 RESULTS 


4.1 Performance 


In this section we test the performance of the emacss-mcmc 
method. We first determined a suitable number of the walk¬ 
ers, n w , burn-in iterations, n b, and subsequent chain itera¬ 
tions, n c , such that we have a good balance between proper 
coverage and quick convergence. To this end we ran a dozen 
simulations for the cluster M4 varying these three numbers 
and plotting the posterior probability of an iteration as a 
function of the iteration number. We divided all the iteration 
in bins of 50 iterations and calculated the minimum, max¬ 
imum and mean probability per bin, see Figure [3] In these 
simulations we used a prior distribution which is uniform in 
mass (similar to the first line of Eq. (|3|), but semi-uniform 
in half-mass radius: the upper limit of the half-mass radius 
is mass-dependent through the dependence on the Jacobi 
radius, rj: 

log(M 0 b 3 ) ^ < log(Mobs) + 3, 

0 < < 0.3rj. (9) 


This mass-dependent upper limit of the initial half-mass 
radius is motivated by the fact that most clusters with 
initial half-mass radii larger than 30% of the Jacobi radius 
will quickly dissolve (Alexander et al. 2014). As we show 
later on in this section, initializing the walkers according 
to this prior does not exclude the parameter space with 
rhm /rj > 0.3, but it does accelerate the convergence. 


In order to judge whether the walkers have converged, 
we look at the instance that both the maximum and the 
mean posterior probability have stabilized, which means 
that its value does not increase or decrease by a significant 
amount, say 100%, and thus only show the variation caused 
by the scatter. The scatter in the mean posterior probability 
is a natural feature of any MCMC sampler, because it is 
important that even though a (local) maximum in posterior 
probability has been found, the walkers continue to explore 
other regions of parameter space, in order to locate possible 
other maxima. 


From the first row of Figure [3j showing the simula¬ 
tion with (n w , rib, n c ) = (20, 100, 100), we see that 
the walkers already start to converge after about 2000 
iterations, marking the end of the burn-in phase. However, 
from the right column of the first row, we see that the 
coverage of the M — rhm plane is still poor, meaning that 
we can observe by eye that the two-dimensional parameter 
space is not well-sampled and/or that the sampled region 
does not cover a large region of that parameter space. In 
the second row of Figure [3] we see that both convergence 
and proper coverage, by eye, seem to have been established 
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Figure 3. MCMC performance in terms of convergence for the simulation of M4 without a correction for mass segregation. The left 
column shows the posterior probability of an iteration as a function the iteration number, averaged over bins of 50 iterations. The solid 
black line shows the mean probability of each bin, the pink area marks the region between the minimum and the maximum probability 
in each bin and the vertical solid cyan, dashed blue and dotted red lines, mark the phases after initialization, the burn-in phase and the 
subsequent chain iterations, respectively. The middle column shows the distribution of the walkers after initialization (our prior; cyan 
crosses), burn-in (blue circles) and the chain iterations (red triangles), respectively. The right column shows the distribution of all the 
sampled initial conditions with a non-zero posterior probability, colour coded by the value of the posterior probability as indicated by 
the color bar. The first row is for the simulation with (n w , rib, n c ) = (20, 100, 100), the second row for the simulation with (n w , rib, n c ) 
= (50, 100, 1000) and the third row for the simulation with (n w , n^, n c ) = (100, 1000, 10000). 

































































The initial conditions of observed star clusters I. 13 


for the simulation with (n w , rib, n c ) = (50, 100, 1000): the 
walkers probed a wide range of initial masses and half-mass 
radii in proposing initial conditions. Nevertheless, we have 
chosen to be conservative and to use (n w , rib, n c ) = (50, 
100, 1000) only for test simulation^ and to use (n w , rib, 
n c ) = (100, 1000, 10000) for all main simulations in this 
work, see the third row of Figure [3] 

We secondly determined whether the choice of a prior 
distribution would effect the ranges of parameter space 
that are covered and thus the determined probable initial 
condition distribution. To test this we ran three simulations 
to determine the initial conditions for the cluster M4 with 
(n w , rib, n c ) = (50, 100, 1000) using three different priors, 
but otherwise being the same. These three priors were: 

1 a uniform distribution in both (logaritmic) mass and 
half-mass radius, according to Eq. 

2 a normal distribution in both parameters with the mean 
equal to log(M 0 bs) and rhm,obs respectively, and a standard 
deviation equal to 10% of the mean; 

3 a semi-uniform distribution, according to Eq. 

The results of these three simulations are shown in Figure [4] 
From this figure we can see that the covered area in the 
M — fhm plane of the probable initial conditions is similar. 
If we consider all sampled initial conditions (also those 
outside of the boundaries given by Eq. e.g. negative 
half-mass radii, with posterior probabilities equal to zero, 
which is not shown in Figure [4|, the sampled area is 
significantly different: the simulation with a uniform prior 
covers a larger range in initial half-mass radius (—500 < 
(rhm,i/pc) < 1000) and slightly different, but overlapping 
range in logarithmic mass (2.0 < log(Mi/M 0 ) < 10.5) com¬ 
pared to the semi-uniform prior (—68.9 < (rh m ,i /pc) < 150 
and 2.9 < log(Mi/M 0 ) < 9.8) and the normal prior ( 
-26.9 < (r hm ,i/pc) < 60.4 and 1.1 < log(Mi/M 0 ) < 9.1). 
If we consider initial conditions with posterior probability 
greater than zero, then the covered area is slightly more 
similar: the walkers initialized according to a uniform prior 
reach larger initial half-mass radii, but the initial mass 
range is comparable (uniform: 10" 3 < (rh m, i /pc) < 95.7 
and 5.3 < log(Mi/M 0 ) < 7.8; semi-uniform: 10 -4 < 
(rhm,i/pc) < 63.8 and 5.3 < log(Mi/M 0 ) < 7.8; normal: 
10 -4 < (rhm,i/pc) < 58.0 and 5.3 < log(Mi/M 0 ) < 7.7). We 
observe in the top panel of Figure [4] that there is an area 
of initial conditions with high initial mass and high initial 
radius, which is not probed with the simulations with the 
other two priors. However, this is not a favourable region 
in terms of posterior probability and thus we conclude 
that simulations with different prior distributions properly 
cover the relevant ranges of parameter space and that the 
derived initial conditions are prior independent. Similar 
behaviour is seen for the other clusters in our sample. We 
use the semi-uniform prior distribution for the rest of our 
simulations in this work. 

The final performance characteristic to test is the sampling 
of initial conditions by the walkers. The results of all the 


6 Testing the effect of the prior distribution or observational er¬ 
rors on the determined probable initial condition distributions. 


simulations for both the 2D and 5D runs are shown in 
Figures |6][T4| Figure [T5] and Figures |AT||A8| By comparing 
the two-dimensional histograms with the confidence con¬ 
tours in these figures, we see that the most sampled areas 
overlap with the high posterior probability regions. This 
is what one would expect for an MCMC method with a 
sufficient number of iterations, and hence shows the proper 
performance of our method. We note that the number of 
sets of initial conditions in the 99.7% and 68.3% confidence 
regions is less than ~99.7%, respectively ~68.3%, of the 
total number of surviving initial conditions, and that 
increasing the number of iterations does not increase these 
percentages. For instance, for M4 the number of sets of 
initial conditions in the 99.7% and 68.3% confidence regions 
is ~80% respectively ~40% of the total number of surviving 
initial conditions for both the simulation with (n w , n b, n c ) 
= (50, 100, 1000) and the simulation with (n w , rib, n c ) = 
(100, 1000, 10000). This is again because even though the 
walkers have converged to high posterior probability regions 
of parameter space, they continue to sample unexplored 
regions as well. 


4.2 Direct comparison 

Figure [5] shows the direct comparison between emacss and 
the dedicated studies by running emacss from their best-fit 
initial condition. All results presented in this section are 
summarized in Table [3j where we compare the best-fit 
results of the dedicated study (DS) and emacss’ result 
when starting from the initial condition of the dedicated 
study in row two and three for each cluster. 

From Figure [5] we see that the emacss results com¬ 
pare quite well to the Monte Carlo and the iV-body results 
for NGC 6397, M4, M22, Pal 4, 47 Tuc, Pal 14 and Gl, 
with a difference in final conditions < 25 % with respect to 
the dedicated studies in both (linear) mass and half-mass 
radius for each of these clusters. For M4, M22 and 47 
Tuc, we see that the evolution with emacss led to a 
close match in final conditions (with < 8% difference in 
both mass and half-mass radius) losing slightly more mass 
during its evolution, and reaching smaller final radii. For 
these three clusters the final radii are smaller, because 
the cluster did not expand as much as in the MOCCA 
simulation. For M4 the radius did expand up to 2.87 pc 
in the evolution with emacss, but near the end of the 
simulation the cluster already started to contract due to 
stellar evaporation. For NGC 6397 we have almost an 
exact match in final half-mass radius (< 1% difference), 
but the amount of mass loss is less, leading to a difference 
of ~ 21 % in linear mass with respect to the dedicated study. 


For Gl the cluster modeled with emacss more mass 
than the cluster modeled with a scaled iV-body simulation 
(~ 6% difference in mass), but it also expanded more 
(~ 25% difference in radius). However, we note that in 
calculating the scaled-up version of the initial mass and 
half-mass radius of Gl found by Baumgardt et al. (2003b), 
we assumed that the same amount of mass loss and stellar 
evolution induced expansion occured as in the small scale 
model. This does not have to be the case. Therefore it could 
be that the scaled up initial condition from which we start 
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Figure 4. Similar to Figure[3] but this time showing the dependence of the initial conditions and performance in terms of convergence and 
coverage on the choice of prior distribution. The top panel shows the performance after choosing a uniform prior distribution indicated by 
the square box. The middle panel shows the performance after choosing a normal prior distribution, with the current mass and half-mass 
radius as the mean and a tenth of these values as the standard deviation; see the ellipses indicating the one sigma (solid black line) and 
the three sigma lines of the prior distribution in the zoomed in smaller panel. The third panel shows the performance after choosing a 
semi-uniform distribution, i.e. the mass is drawn uniformly between two boundaries, but the radius is dependent on the Jacobi radius, 
which is mass dependent; see the semi-uniform prior indicated by the black dashed line. 
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Figure 5. Direct comparison between emacss and the dedicated studies (DSs) by running emacss from the best-fit initial condition 
of each DS. Each panel shows the total mass versus half-mass radius for one of the nine validation clusters. The blue diamond(s) and 
cross(es) show the best-fit initial and final condition(s) respectively of the DS; the black dashed line shows the evolution with emacss 
when started from the initial condition of the DS and the black square shows the final condition of this evolution with emacss. For uj Cen 
we show a blue line, indicating the final mass of the DS, since their final half-mass radius was not given. For Pal 4 and Pal 14 we compare 
to more than one model, so we connect the final conditions of the DS and emacss with a red solid line for clarity. As a comparison we 
also plot th e minimum densities a cluster needs to have in order to be stable against the tidal disruption of a galax y (phm ~ 0.1 M©pc -3 
(Bok 1934|, green triangle line) and against passing giant molecular clouds (phm ~ 10 M©pc -3 ( jSpitzer 1 958| ), cyan star line). The 
turquoise hexagon line shows a mean density of ~ 5 • 10 7 M©pc -3 that might be required for a cluster with a half-mass radius of 0.2 pc 
to form an IMBH via a runaway merger ( |McMill an 2008| >. We furthermore show the range of observed half-mass densities for globular 
clusters, phm ~ 10 -2-3,7 M©pc -3 , ( red shaded regio n ) and the range of present-day observed half-mass densities for clusters younger 
than 10 Myr (phm ~ 10 2-4 Mqpc - 3 (Portegies Zwart et al. 2010) (blue shaded region); the purple shaded region is the overlap of these 
blue and red regions. 
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evolving with EMACSS is different from the actual scaled 
up version of the initial condition of Baumgar dt et al.| 
(2003b), which would cause the differences in final mass 


and half-mass radius. For Pal 4 and Pal 14, we see that on 
average the direct iV-body modeling caused the modeled 
clusters to lose slightly more mass than the clusters modeled 
with emacss leading to < 15% difference in mass for Pal 4 
and < 21 % difference in mass for Pal 14. For both Pal 14 
and Pal 4, the emacss clusters expanded a bit more (< 4% 
and <2% difference in radius, respectively). 

For uj Cen the emacss results compare less well to 
those obtained from MOCCA. Their simulated cluster lost 
substantially more mass (leading to a ~ 59 % difference 
in linear mass). We could not compare the final half-mass 
radii, since Giersz & Heggie (2003) did not provide this, as 
described in Section |3.2.1.4| It could be that our assumption 
for calculating what the initial half-mass radius for the 
best-fit of Giersz & Heggie ( 2003| (see Section 3.2.1.4) 
was incorrect and that their initial half-mass radius was 
somewhat smaller or larger. Starting at a different initial 
radius could lead to a different amount of mass loss. 
However, we have done several emacss runs for uj Cen, 
starting at the same initial mass, but with different initial 
half-mass radii in the range 0.01 - 30 pc and we saw that 
the amount of mass loss changed only minimally. The mass 
loss changed at most by a factor ~ 1.08 between two radii 
in this range. For initial half-mass radii > 20 pc the amount 
of mass loss decreased. Only for clusters initially more 
compact than 0.5 pc will the amount mass loss increase 
significantly, but not as much as in the MOCCA simulation. 
It is also not likely that u Cen started out so compact, see 
Section [473] There is also poor agreement between emacss 
and the direct iV-body results for the open cluster M67. 
Again, we see that the cluster simulated with emacss lost 
substantially less mass. Furthermore, the emacss cluster 
expanded by almost a factor two during its evolution, 
whereas the final half-mass radius of the cluster modeled 
with direct iV-body integration is even slightly smaller than 
the initial one. These differences require some explanation. 


Besides starting from the same initial total mass and 
half-mass radius, we kept the conditions of our simulation 
as similar as possible to those of the dedicated studies. 
However, there are a number of parameters that could not 
be taken equal between the codes. One of these parameters 
is the initial number of stars. Since emacss is currently 
tested against 7V-body simulations for only one value of 
the initial mean mass (m = 0.64 M©) we always used this 
value for all our simulations. This means that once we 
set the initial total mass of the cluster, we immediately 
set the initial total number of stars as well, which was 
different from the initial number of stars in each of the 
dedicated models, see Table [3] This lead to different initial 
half-mass relaxation time scales, t r h, with a larger t r h 
for larger N. The initial half-mass relaxation time in the 
simulations of the dedicated studies for M67, uj Cen, Pal 4 
and Pal 14 respectively was a factor 1.23, 1.47, 1.28 and 
1.12 respectively smaller than the EMACSS simulations. 
One would thus intuitively expect the clusters with shorter 
initial half-mass relaxation times to dissolve quicker, and 
this would lead to a relatively quicker mass loss in the first 


phase of its life ( Lamers et al.p OlO). Furthermore, Heggie 


& Hut (2003) explain that in larger models, i.e. models 


with larger N , the escape rate per relaxation time is larger. 
This could explain why the dedicated studies on M67, uj 
Cen, Pal 4 and Pal 14 lost more mass, since they start with 
a larger number of stars initially. For the clusters M4, M22 
and 47Tuc the evolution with EMACSS started out with 
larger number of stars, and thus larger initial half-mass 
relaxation times by a factor 1.09, 1.07 and 1.25 respectively. 
Here we thus see that the simulated clusters with EMACSS 
lost more mass, but not by much, see Figure [5] For NGC 
6397 and G1 the initial number of stars of the dedicated 
studies was not given. 


Also, the prescription of stellar evolution and the (mass 
limits of the) initial mass function used by emacss and by 
the dedicated studies and sometimes the Galactic potential, 
setting the tidal field, were different. All this taken together 
led only to moderate differences in the final cluster total 
mass and half-mass radius for the clusters NGC 6397, M4, 
M22, Pal 4, 47 Tuc, Pal 14 and Gl, but to more significant 
difference for uj Cen and M67. For M67, for example, it led 
to a larger initial Jacobi radius (rj = 37.6 pc) in EMACSS 
than the (similar) tidal radius (r*t = 31.8 pc) in the N-body 
simulation^ This, in part, explains the smaller amount of 
mass loss in the simulation with EMACSS and hence the 
overall larger amount of expansion. However, for uj Cen 
the difference in mass loss cannot be explained by this, 
since we had an initial Jacobi radius similar to the tidal 
radius in the MOCCA simulation (rj = 89.9 pc; r t = 90 pc). 


Another difference is that all dedicated studies (except for 


Giersz & Heggie 2003) include some direct prescription 


for binaries, whereas the version of emacss used for this 
study does not. For each globular cluster in our sample, the 
dedicated studies had only a small (< 1%) or no primordial 
binary fraction, so the effects that binaries have on the 
evolution of global cluster parameters such as total mass 
and half-mass radius are expected to be small here. Thus 
for validation purposes, binaries are not expected to cause 
major differences in the results between emacss and the 
dedicated studies for this sample of globular clusters. For 
open cluster M67, however, a large initial binary fraction of 
50% was assumed in Hurley et al. (2005) and they showed 
that this fraction even increased throughout the 4 Gyr of 
evolution, also due to the evaporation of single stars. So for 
M67 we do expect binaries to play an important role in at 
least the cluster dynamics and this could certainly result in 
a very different evolution when compared to the evolution 
without binaries. This therefore also contributes to the dif¬ 
ference between emacss and the direct iV-body simulation 
of M67. Just like mass loss due to stellar evolution, hard 
binaries are expected to cause expansion of the half-mass 
radius (Giersz & Heggie 2011), and so based on this type of 
binaries alone one would not expect the cluster evolved with 
emacss to have expanded more. However, it is not obvious 
what to expect in terms of expansion of the half-mass radius 
for a steady 50% fraction of primordial and dynamically 
formed, hard and soft binaries. An interesting side note is 


7 Note again that rj = rt for the type of potential we use here. 
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that when we continued the evolution for M67 for another 
~ 5 Gyr, it reaches a mass and half-mass radius within 4% 


of those of Hurley et al. (2005). This might suggest that the 
presence of binaries accelerates the evolution. To test this 
hypothesis, we are running an N-body simulation similar 


to the one conducted by Hurley et al. (2005), but without 


binaries, of which we will show the results in forthcoming 
work. 


The fact that emacss is designed for evolving globu¬ 
lar clusters instead of open clusters is not in itself a reason 
for the sparse agreement. The majority of open clusters 
reside in the disk and they are more likely to undergo en¬ 
counters with giant molecular clouds, which have dramatic 


effects on the cluster evolution (Gieles et al. 2006). How¬ 


ever, even though emacss does not include these effects, 


neither did the simulations of Hurley et al. (2005), so for 


the comparison between the two codes for the validation, 
this should not matter. Applying our method to an open 
cluster for actually constraining its initial conditions, this 
will become important and some prescription taking into 
account these disruptive events is required. 


4.3 Independent emacss-mcmc runs 


The results of the independent emacss-mcmc runs in 2D 
are shown in Figures [6] to |14| and are summarized in Table [3j 
where we compare the observations and emacss’ best-fit re¬ 
sults for the simulations where we did not correct the obser¬ 
vations for mass segregation (see Section 2.3) and emacss’ 
best-fit results for the simulations where we did correct for 
mass segregation, in row one, four and five, respectively, per 
cluster. We see that for some clusters (M67, NGC 6397, M4, 
M22 and 47 Tuc) the initial conditions extend to small half¬ 
mass radii and high initial half-mass densities. These initial 
conditions correspond to an average distance between the 
cluster stars on the order of 10 AU, and thus these densities 
are too high to be physical. In Figures [6] to [T4] we show the 
range of observed half-mass densities for globular clusters, 
phm ~ 10 -2-3 ' 7 M 0 pc -3 , calculated according to 


/?hm = 3 Mj (8-7rrh m ), (10) 



shaded region). The high-phm initial conditions mentioned 
above are significantly larger than the ones observed today 
and this is attributable to the small initial half-mass radius 
the MCMC sampled for these clusters. One could argue that 
the current half-mass densities need not be representative 
for the initial half-mass densities, and (globular) clusters 
may have had larger initial half-mass densities at younger 
ages. However, from Figures [6] to |9] we can see that the 
range of present-day observed half-mass densities for clus¬ 
ters younger than 10 Myr is still a few orders of magnitudes 
less dense than the initial densities for M67, NGC 6397, M4 
and M22 in our calculations. One could also argue that the 
precursors of the old (globular) clusters might have been 


very different from the young clusters today. High densities 
might even be essential to enable runaway mergers as a 
pathway to produce intermediate black holes (IMBHs) in 
globular clusters (McMillan 2008). McMillan (2008) gives 
the example that a mean density of ~ 5 • 10 Y M 0 pc would 
be required for a cluster with a half-mass radius of 0.2 pc. 


Moreover, Pfalzner et al. (2014) ran semi-analytical models 


for the formation of star clusters and show that the central 
cluster area can have stellar densities of ~ 4-10 5 M 0 pc -3 at 
t he m oment of gas expulsion; see Figure 2 of | Pfalzner et al.] 
(2014). However, densities greater than ~ 10 10 M 0 pc -3 
seem too extreme. Moreover, it is important to keep in 
mind that the initial conditions we derive in this work are 
the conditions of a star cluster after residual gas expulsion 
and re-virialisation. Since all star clusters expand due to 


residual gas expulsion (Baumgardt & Kroupa 2007), the 


clusters are expected to be even denser, i.e. more massive 
and more compact, directly after cluster formation. 


The reason the MCMC code selected these initial conditions 
is that we have not build in a criterium for only selecting 
initial conditions below a maximum allowed initial half-mass 
density, simply because we do not know what this upper 
limit should be. We also considered it to be better not 
to limit the MCMC code, but to include possible density 
limits only in the analysis phase. Furthermore, for the four 
clusters (M67, NGC 6397, M4 and M22) where probable 
initial conditions are found in the high half-mass density 
regions of parameter space, we see that the probable initial 


conditions are degenerate (see Section 4.3.3). From the 
brown contours in Figures [6] to [9] we can conclude that lower 
density initial conditions are practically equally probable. 


The results of the 5D simulations are shown in Fig¬ 
ure [15] and in the figures in the Appendix [A] Table [4] 
furthermore shows the characteristics of our best-fit model 
without a correction for mass segregation. From the brown 
contours in Figure [15] we can once again see that the most 
probable 0.1% of the initial conditions are found in the 
most sampled region, but that a probable initial condition 
does not have to have each parameter to originate from the 
most sampled region. 


4-3.1 Correcting observations for mass segregation 


The dedicated studies to which we compare our results 
did not correct their observations for mass segregation. 
This is because it is difficult to determine the amount of 
mass segregation that the cluster experienced and hence 
the ratio rhm/'Tphi at f = r Q bs, i.e. the factor cms in 
Section [T3] However, many globular clusters show signs of 
mass segregation. Even the outer halo cluster Pal 14 has 
recently been shown to be (primordially) mass segregated 
(Frank et al. 2014). For this reason we investigated what 


the influence of observing a mass segregated cluster, but 
not correcting the half-mass radius for it, could have on the 
determined initial conditions of the cluster. We did this for 
a fixed ratio 7*hmA*phi, see Section : 


2.3 


but the outcome will 
illustrate the importance of correcting for mass segregation 
on the interpretation of the evolution of the observed cluster. 


The first thing we can see from Figures [6] - [14] is that 
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Figure 6. The same as in Figure[2] but for the cluster M67 and with more details: the contours and histograms of the 99.7% confidence 
level (dotted), the 68.3% (solid) confidence level and the p > 0.9p ma x (solid) initial conditions for both the simulations without (yellow, 
red and brown, respectively) and with (cyan, blue and dark blue, respectively) a correction for mass segregation. The green contour and 
histograms show the clusters which dissolved before reaching the age of the cluster, r Q b s - The shaded regions, the lines depicted by green 
triangles, cyan stars and turquoise hexagons are as in Figure [5] 


correcting the observations for mass segregation or not, 
affects the shape and location of the corresponding 99.7% 
confidence region in the initial conditions. In other words, 
it changes which initial conditions are considered probable. 
In some cases these changes are mild, such that the 99.7% 
confidence regions still show some overlap, e.g. for M67, 
NGC 6397, M4, M22; in Section |4.3.3| we will see that 
this has to do with the degeneracy due to core-collapse. 
But in the other cases the shapes or location of the 99.7% 
confidence regions are different. And in almost all cases, the 
contours of highest probability are completely separated. 
It even turns out that the factor cms = 1.9 we used for 
the globular clusters (Hur ley|[2 007) is typical for clusters of 
low density, whereas this factor seems to increase for more 
dense clusters (Mirek Giersz, private communication). For 
denser clusters we thus expect the differences in the initial 
conditions for the corrected and non-corrected simulations 
to be even more distinct. However, many studies do not aim 
to (just) fit structural parameters such as half-mass radius 
and total cluster mass, but aim to (also) have a good fit in 
the surface brightness profile, velocity dispersion profile and 
luminosity profile at multiple radii (|Giersz &; Heggie|[ 


Hurley et al.|2005||Heggie &"( 4iersz| 


Baumgardt et al. 


2008[ Giersz & Heggie 2009[ |2011[ |Heggie &; Giersz||2014| ); 

see Section |4.4| It wilTtherefore be interesting to investigate 

whether a cluster that will be modeled with iV-body or 
Monte Carlo simulations starting from one of the best-fit 
initial conditions with a correction for mass segregation, 
can produce the three above-mentioned profiles that match 
the observed ones equally well as a cluster starting at the 


best-fit initial conditions without this correction. 

For further validation purposes, we focus on the com¬ 
parison between our simulations without a correction for 
mass segregation to the results of the dedicated studies. 


4.3.2 Dissolving clusters 

From the Figures [6] to [14] we see that for each of the nine 
validation clusters, the modeled clusters which do not sur¬ 
vive up to the cluster age are in most cases relatively larger 
(and sparser) and/or less massive than the clusters which do 
survive. Sometimes these clusters dissolved too sparse initial 
conditions were chosen by the MCMC code, which led to im¬ 
mediate (within a few Myr) dissolution. E.g. this is the case 
for all dissolving clusters for Pal 14, which have phm >0.1 
M oP c- 3 . In other cases these clusters underwent relatively 
quick mass loss, dissolving before reaching the observed clus¬ 
ter’s age. The surviving and dissolving initial conditions are 
well separated, except at the borders of the two distribu¬ 
tions, which is to be expected. From Figure [T5| we see that 
for the 5D runs the dissolution regions projected onto a two- 
dimensional plane are not well separated from the surviving 
(high posterior probability) regions anymore. See e.g. panel 
11 of Figure [l5] where the 99.7% confidence contour largely 
overlaps with the dissociation region contour. The reason 
for this overlap is simply because that figure is a projec¬ 
tion of a five-dimensional space of initial conditions onto 
a two-dimensional space. In other words, a suitable initial 
condition (rh m ,i, Mi) for one combination of the nuisance pa- 
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Figure 7. The same as Figure [6] but for the cluster NGC6397. 



Figure 8. The same as Figure [6] but for the cluster M4. 


rameters ( Rgc , v, r), could be dissolving for another combi¬ 
nation of these nuisance parameters. If we were to take 2D 
projections of this five-dimensional space at distinct combi¬ 
nations of (.Rgc,^, t), the dissociation regions should again 
be well separated from the surviving regions, except at bor¬ 
ders where we expect to have an overlap again. 


4.3.3 Degeneracies: dynamical age and morphology 

We can divide the surviving clusters in two types of clusters: 
1) the clusters which underwent a core-collapse and are 
now in the third and ‘balanced’ phase of their evolution 
and 2) the clusters which did not undergo a core-collapse 
(yet) and are still in the second and ‘unbalanced’ phase 
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Figure 9. The same as Figure [6] but for the cluster M22. 



Figure 10. The same as Figure [6] but for the cluster Pal4. 


of their evolution (see Section 2.1 for a definition of these 
evolutionary phases). See e.g. Figure 16 where we make 
a distinction between the initial and final conditions of 
the model clusters which undergo a core-collapse during 
their evolution (red) and those which do not (blue) for the 
2D run for the cluster M4 without a correction for mass 


segregation. We can see that the clusters which undergo 
a core-collapse in our simulations generally start out with 
smaller half-mass relaxation times, log(t r h ) < 2.8 in this 
example for M4; these clusters also have relatively small 
radii. For those clusters with log(t r h) > 2.8 we see that 
for a fixed radius the clusters that did not yet undergo a 
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Figure 11. The same as Figure [6] but for the cluster 47Tuc. 



Figure 12. The same as Figure [6] but for the cluster Pall4. 


core-collapse are more massive. This is what one would 
expect from the fact that the core-collapse time, t cc , scales 
with the half-mass relaxation time, t r h, which increases 
for large rhm and large N. The lower mass clusters with 
log(t r h) > 2.8 also form a degeneracy towards larger radii 
and masses. With EMACSS we keep track of the evolution 


of the derivative of the half-mass radius, drhm/dt. It turns 
out that the degeneracy towards large radii is produced 
by the fact that those clusters experienced a phase of 
contraction. This is consistent with the fact that once 
clusters have reached a certain relatively large radius (due 
to expansion), the escape of stars over the Jacobi radius 

































































log(M/M 0 ) log(M/M 


22 J. T. Pijloo et al 



Figure 13. The same as Figure [6] but for the cluster uo Cen. 



Figure 14. The same as Figure [6] but for the cluster Gl. 


becomes dominant (Gieles & Baumgardt 2008). The large 
amount of mass loss of these clusters corresponds to the 
decrease of their half-mass relaxation time, such that they 
will still undergo a core-collapse. This degeneracy is limited 
on the low density side by the dissolution region. 


When we examine Figures [6] to [T4] we notice that the 
99.7% confidence regions in the initial conditions are built- 
up out of three characteristic shapes (see also Figure [l7|): 

(i) a roughly horizontal distribution with a negative slope 
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Cluster 

Source 

Ni 

Mi 

[10 5 M 0 ] 

^ hm,i 

[pc] 

N / 

M/ 

[10 4 M@] 

r/im,/ 

[pc] 

CC? a 

M67 

Observation 

- 

- 

- 

- 

0.140 

3.35 / 4.52 

- 


dedicated study 

36.0k 

0.187 

3.90 

3.52 k 

0.204 

3.80 

no 


EMACSS from IC of DS 

29.2k 

0.187 

3.90 

19.0 

0.820 

6.51 

no 


EMACSS best-fit 

25.2 k 

0.162 

0.006 

2.64 

0.141 

3.39 

yes 


EMACSS best-fit (MS corr.) 

347 k 

2.22 

43.0 

2.39 

0.149 

4.49 

no 

NGC 6397 

Observation 

- 

- 

- 

- 

6.60 

2.31 / 3.29 

yes 


dedicated study 

- 

3.65 

0.400 

- 

6.03 

3.22 

yes 


EMACSS from IC of DS 

570 k 

3.65 

0.400 

143 k 

7.28 

3.21 

yes 


EMACSS best-fit 

1.21 M 

7.77 

1.18 TO -4 

102 k 

4.40 

2.89 

yes 


EMACSS best-fit (MS corr.) 

435 k 

2.79 

2.78 

138 k 

6.60 

3.29 

yes 

M4 

Observation 

- 

- 

- 

- 

6.30 

2.44 / 3.48 

no 


dedicated study 

485 k 

3.40 

0.580 

86.0k 

4.61 

2.89 

yes 


EMACSS from IC of DS 

531k 

3.40 

0.580 

79.1k 

4.54 

2.69 

yes 


EMACSS best-fit 

1.20 M 

7.66 

1.19 TO -3 

121k 

5.87 

2.78 

yes 


EMACSS best-fit (MS corr.) 

495 M 

3.17 

4.08 

125 k 

6.28 

3.48 

no 

M22 

Observation 

- 

- 

- 

- 

33.0 

4.17 / 5.95 

no 


dedicated study 

832 k 

5.70 

2.72 

- 

32.0 

6.60 

no 


EMACSS from IC of DS 

891k 

5.70 

2.72 

811k 

29.9 

6.28 

no 


EMACSS best-fit 

1.17M 

7.50 

0.358 

849 k 

33.1 

4.19 

yes 


EMACSS best-fit (MS corr.) 

978 k 

6.26 

2.53 

898 k 

33.1 

5.95 

no 

Pal 4 

Observation 

- 

- 

- 

- 

2.98 

24.5 / 35.0 

no 


dedicated study 

100 k 

0.500-0.570 

12.0-14.5 

- 

2.68-3.24 

23.2-27.6 

no-no 


EMACSS from IC of DS 

78.1-87.3 k 

0.500-0.570 

12.0-14.5 

76.5-87.2 k 

2.78-3.18 

26.6-31.1 

no 


EMACSS best-fit 

83.7k 

0.535 

10.9 

82.1k 

2.98 

24.5 

no 


EMACSS best-fit (MS corr.) 

83.8 k 

0.536 

16.6 

81.9k 

2.98 

35.0 

no 

47 Tuc 

Observation 

- 

- 

- 

- 

110 

4.87 / 6.94 

no 


dedicated study 

2.04 M 

16.4 

1.91 

1.87M 

90.0 

4.96 

no 


EMACSS from IC of DS 

2.56 M 

16.4 

1.91 

2.38 M 

87.4 

4.57 

no 


EMACSS best-fit 

3.19 M 

20.4 

2.04 

3.01 M 

110 

4.87 

no 


EMACSS best-fit (MS corr.) 

3.22 M 

20.6 

3.12 

3.00 M 

110 

6.93 

no 

Pal 14 

Observation 

- 

- 

- 

- 

1.1725 

35.2 / 50.2 

no 


dedicated study 

70.0-100 k 

0.400-0.600 

15.0-25.0 

- 

1.86-2.89 

26.1-42.8 

- 


EMACSS from IC of DS 

62.5-93.8 k 

0.400-0.600 

15.0-25.0 

59.4-87.9 k 

2.17-3.22 

31.5-47.9 

no - no 


EMACSS best-fit 

35.6k 

0.228 

17.6 

31.5k 

1.17 

35.2 

no 


EMACSS best-fit (MS corr.) 

136 k 

0.869 

64.1 

56.2 k 

2.11 

49.3 

no 

lj Cen 

Observation 

- 

- 

- 

- 

390 

9.89 / 14.1 

no 


dedicated study 

- 

110 

9.33 

- 

360 

- 

no 


EMACSS from IC of DS 

17.2 M 

110 

9.33 

15.7M 

570 

16.5 

no 


EMACSS best-fit 

11.7M 

74.7 

5.14 

10.6 M 

390 

9.89 

no 


EMACSS best-fit (MS corr.) 

11.9M 

76.2 

7.84 

10.7M 

390 

14.1 

no 

Gl 

Observation 

- 

- 

- 

- 

1.50 TO 3 

17.6 / 25.1 

- 


dedicated study 

- 

127 

3.86 

- 

760 

6.76 

no 


EMACSS from IC of DS 

19.8 M 

127 

3.86 

19.8 M 

713 

8.5 

no 


EMACSS best-fit 

41.7M 

267 

9.20 

41.7M 

1.50 TO 3 

17.6 

no 


EMACSS best-fit (MS corr.) 

41.7M 

267 

13.5 

41.7M 

1.50 TO 3 

25.1 

no 


Table 3. Results for the validation in three significant figures. For each validation cluster mentioned in the first column seven parameters 
are compared between the observation in the first sub-row, the best-fit results of the dedicated study (DS, see Table [l] for the references) 
in the second sub-row, emacss’ result from the initial condition of the DS in the third sub-row, emacss-mcmc’ best-fit results without 
and with a correction for mass segregation in sub-rows four and five respectively (see Section |2.3| . The parameters are the initial number 
of stars Ni in column three, the initial mass (in 10 5 M®) in column four, the initial half-mass radius rh m ,i (in pc) in column five, the 
final number of stars Nf in column six, the final (current) mass M f (in 10 4 M®) in column seven, the final (current) half-mass radius 
fhm,f (in pc) in column eight and the question whether mass segregation has occurred in column nine. The observed value for the current 
half-mass radius in column eight shows two values: the first one without a correction for mass segregation and the second one with a 
correction for mass segregation, see Section [2. 3 1 and Table [2] 

a The results whether a cluster had undergone core-collapse according to observations are taken from Trager et al.| (1995]>, except for 
Gl, where we adopt what [Baumgardt et al.] ( |2003b| > argued about the core-collapse state of these two clusters. 
b The authors studying this cluster with their dedicated study find multiple initial conditions; we only list their first model here. 


towards small half-mass radii. This distribution has high 
posterior probabilities for a large spread in half-mass radii 
and a moderate spread in mass. The 99.7% confidence re¬ 
gions of the clusters M67, NGC 6397, M4 and M22 have 
this morphology for both our models with and without a 


correction for mass segregation (MS). This distribution cor¬ 
responds to the clusters that undergo a core-collapse within 
their evolution time, r 0 bs- 

(ii) a distribution with a positive slope reaching to large 
radii. This characteristic shape is seen in the initial condi- 
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Figure 15. Initial (and final) condition distributions for cluster M67 from the 5D MCMC simulations with fitting parameters mass and 
half-mass radius, and nuisance parameters Galactocentric radius, orbital velocity and age, each with 10% errors for the simulations which 
fit the observables without a correction for mass segregation. Numbered from top to bottom, from left to right, in panels two, four, five, 
seven, eight, nine, eleven, twelve, thirteen and fourteen show two-dimensional histograms in a black-white density plot of two of the five 
parameters initial mass, initial half-mass radius, initial Galactocentric radius, initial orbital velocity and (final) age against one another: 
the darker the area, the more initial conditions in this area were sampled. Over-plotted on these two-dimensional histograms are the 
outer contours of the 99.7% (yellow) and 68.3% (red) confidence contours for the simulation fitted to observations without a correction 
for mass segregation. The cyan and blue contour in panel eleven show the 99.7% and 68.3% confidence contours of the final conditions in 
total mass vs half-mass radius, respectively. The brown contours show the most probable region with p > 0.999p m ax- The green contours 
show the clusters which dissolve before reaching the age of the cluster, r Q b s - The red square with error bars show the observation of the 
cluster’s current mass and half-mass radius when no correction for mass segregation is made, and the blue square with error bars shows 
the observable with a correction. The pink filled circle denotes the initial condition used by the dedicated study, which evolves to the 
final condition shown by the grey filled circle in the right panel. The corresponding projected histograms are shown in panels one, three, 
six, ten and fifteen. 


tions for the clusters M67, NGC 6397, M4 and Pal 14 for 
the models with and without a correction for MS and for uj 
Cen for the models with a correction for MS. This distribu¬ 
tion corresponds to the clusters that experienced a phase of 
contraction. 

(iii) a two-dimensional gaussian-shaped distribution at 
relatively larger initial half-mass radii than the previous 
case. The 99.7% confidence regions of the clusters Pal 4 and 


G1 have purely this morphology for the models with and 
without a correction for MS. This shape is also seen in the 
99.7% confidence regions of the clusters M22, 47 Tuc, Pal 14 
and uj Cen for the models with and without a correction for 
MS. The initial conditions from this distribution correspond 
to clusters that do not undergo a core-collapse within the 
evolution time. In this case the initial conditions are better 
constrained. 
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Cluster 

Ni 

Mi 

[10 5 M 0 ] 

^ hm,i 

[pc] 

Prg 

kpc 

V 

km/s 

N / 

M f 

[10 4 M 0 ] 

r hm,f 

[pc] 

t 

Gyr 

M67 

17.0k 

0.108 

1.24 

7.49 

216 

2.18k 

0.139k 

3.35 

4.93 

NGC 6397 

579 k 

3.70 

1.52 

1.51 

265 

99.5 k 

6.37 

2.37 

9.89 

M4 

820 k 

5.25 

0.179 

1.45 

244 

98.3 k 

6.29 

2.45 

11.5 

M22 

988 k 

6.32 

1.65 

4.59 

197 

516 k 

33.0 

4.17 

10.8 

Pal 4 

83.4 k 

0.534 

11.1 

96.6 

219 

46.6 k 

2.98 

24.5 

10.6 

47 Tuc 

3.21 M 

20.5 

1.97 

3.57 

219 

1.72 M 

110 

4.87 

13.8 

Pal 14 

34.8 k 

0.222 

16.8 

88.3 

236 

18.3k 

1.17 

35.2 

12.1 

uj Cen 

11.7M 

74.9 

5.12 

1.23 

218 

6.09 M 

390 

8.89 

12.4 

G1 

41.2M 

264 

9.38 

42.8 

234 

23.4 M 

1501 

17.6 

10.3 


Table 4. Best-fit parameters for the 5D simulations fitting the observations without a correction for mass segregation. For each validation 
cluster mentioned in the first column nine parameters of the best-fit model given. The parameters are the initial number of stars Ni in 
column two, the initial mass (in 10 5 Mq) in column three, the initial half-mass radius rh m ,i (in pc) in column four, the Galactocentric 
radius Rgc (in kpc) in column five, the orbital velocity v (in km/s) in column six, the final number of stars Nf in column seven, the 
final mass M f (in 10 4 Mq) in column eight, the final half-mass radius rh m ,f (inpc) in column nine and the age t of the cluster (in Gyr) 
in column ten. 
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Figure 16. The distinction between the initial conditions (top 
panel) and final conditions (bottom panel) of the model clusters 
which undergo a core-collapse during their evolution (red) and 
those which do not (blue) for the 2D run for the cluster M4 with¬ 
out a correction for mass segregation). 


We see the morphology of 99.7% confidence regions 
are built-up out of one or more of the three shapes above. 
The 99.7% confidence regions of the clusters M67, NGC 
6397 and M4 are built-up out of morphology (i) and (ii) 
for models with and wqithout a correction for MS. The 
99.7% confidence regions of the clusters M22 for both 
MS-corrected and not MS-corrected models and 47 Tuc for 
the not MS-corrected models have a morphology consisting 



Figure 17. A schematic overview of the three morphologies: (i) 
a roughly horizontal distribution with a negative slope towards 
small half-mass radii; (ii) a distribution with a positive slope 
reaching to large radii; (iii) a two-dimensional gaussian-shaped 
distribution at relatively larger initial half-mass radii than the 
previous case. 


of shapes (i) and (iii), which shows that both core-collapsed 
and not core-collapsed clusters could fit the observed mass 
and half-mass radius. The 99.7% confidence region of the 
clusters Pal 14 for both MS-corrected and not MS-corrected 
models and u Cen for models with a correction for MS 
consist of morphological shapes (ii) and (iii). 


From the morphologies (i) and (iii) in Figures [6] to 
|14| we see that the morphology roughly changes along 
the sequence (i) —> (i) & (iii) —> (iii). Since the clusters 
are mentioned in the order of increasing current half-mass 
relaxation time, which is a proxy for the dynamical age of 
a star cluster, we thus see a connection between the mor¬ 
phology of probable initial conditions and the dynamical 
age of a cluster. The connection is quite intuitive. Firstly, 
the morphology is connected to the core-collapse state of 
the cluster ((i) - core-collapsed; (iii) - not core-collapsed). 
When a cluster undergoes core-collapse, the half-mass 
radius increases, because energy is injected into the halo 
(Lynden-Bell & Eggleton 1980 Baumgardt et al. 2002). 
The amount by which the half-mass radius increases de¬ 
pends on many different cluster characteristics, such as the 
binary fraction, as well as on the interactions taking place. 
Two clusters with the same initial conditions, but with a 
different initial distribution of positions and velocities of the 
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stars, will have different pathways to core-collapse and the 


halt thereof, also known as statistical fluctuations (Giersz 


& Heggie 2009). It is therefore practically impossible to 


sharply constrain the initial half-mass radius of the cluster, 
because we do not know the amount by which it increased 
during core-collapse. So there are many different possible 
initial half-mass radii and hence we see a degeneracy in the 
initial distribution of half-mass radii. When on the other 
hand the cluster has not yet undergone core collapse, there 
is only a narrow range of possible initial half-mass radii and 
thus the initial half-mass radii are better constrained. 


Secondly, when the half-mass relaxation time is small, 
so is the core-collapse time, since t cc oc thm- So an observed 
cluster with low thm is more likely to have undergone a 
core-collapse and hence its initial condition distribution 
will have morphology (i). If an observed cluster has a 
higher thm, there is a probability that the cluster comes 
from an initial condition which will undergo core-collapse 
during its evolution up to the observed clusters age, but 
also probability that the cluster comes from an initial 
condition which will not. This cluster will thus have an 
initial condition distribution made up of morphology (i) 
and (iii). The clusters with the largest half-mass relaxation 
times will not have initial conditions distributions with 
morphology (i), but rather (iii). 

Morphology (ii) is associated with the contraction of 
the cluster, or, in other words, a decrease of the rhm- And 
just like with the degeneracy associated with core-collapse, 
it is practically impossible to constrain the half-mass radius 
sharply, because we do not know by which amount the 
cluster contracted. So many different evolutionary paths, 
each with different amounts of contraction, could have 
produced the current observables and therefore we have 
this degeneracy. To test our hypotheses for the connection 
between morphology of the initial condition distribution 
and the dynamical age of the cluster, we need to apply 
our method on a larger sample of clusters, which we are 
doing in forthcoming work. If the connection between 
the morphology and the dynamical age is indeed robust, 
then the initial conditions provided by our method gives 
an independent hint onto the core-collapse state of a cluster. 


We will now compare our findings to the observations. For 
M67 we find that our best-fit model without a correction 
for MS underwent a core-collapse. To the best of our 
knowledge M67 is not classified as a post-core-collapse 
cluster and neither have we come across any other claims of 
open clusters going through core-collapse. However, |Hurley| 
et al. (2005) do mention that for binary-rich clusters there 


is no clear core-collapse, at least not in the way that is 
witnessed in simulations without primordial binaries, where 
high core densities need to be reached in order to form 
binaries (see e.g. Hurl ey et al.|[2004 ). The fact that we find 
a post-core-collapse best-fit cluster for M67, is most likely 
a direct result of not directly taking into account binaries 
within emacss, because it is the high binary fraction which 
causes open clusters to not undergo a core-collapse in a way 
that is found for globulars (Hurley et al. 2005). 


Trager et al. (1995) classified uj Cen, M4, 47 Tuc and 


M22 as King profile clusters: these clusters were interpreted 
as clusters which have not undergone core-collapse yet. 
NCG 6397 was classified as a post-core-collapse cluster. 
Defining the term dynamically old (young) for clusters with 
a large (small) age over half-mass relaxation time ratio, 
Tobs/Gh, the cluster NGC 6397 is dynamically old given its 
current half-mass relaxation time of 0.4 Gyr (Harris 2010), 
see also Table [2] Dynamically old clusters are more likely 
to have already ondergone core-collapse, which is consistent 
with the classification. 47 Tuc and uj Cen are dynamically 
young given their current half-mass relaxation times of 3.55 
Gyr and 12.3 Gyr respectively ( Harris| 2010), and are thus 
both far from core-collapse (Trager et al. 1995). However, 


for M22 and M4 Harris (2010) lists relatively short half-mass 


relaxation times, 1.7 Gyr and 0.85 Gyr respectively, such 
that these clusters have undergone a number of relaxation 


times. Moreover, Heggie & Giersz (2008) for the first time 


claimed that M4 is a post-core-collapse cluster, based on 
the observed behavior of the core radius, r c , of their model 
of M4 (see their Figure 5). Therefore, it is interesting to 
see that our best-fit models without a correction for mass 
segregation for both M4 and M22 have also undergone a 
core-collapse, and, as expected, the best-fit model for uj 
Cen did not. 47 Tuc is an interesting case, because we find 
that there are initial conditions of model clusters within 
the 97% confidence region which underwent a core-collapse, 
and initial conditions of model clusters which did not, see 
Section [4. 3.3 1 However, our models show a favor for the not 
core-collapsed solutions, consistent with the observations 
(Trager et al. 1995). 


Both Palomar 4 and 14 are classified as King profile 
clusters, so they are interpreted as clusters which have 


not yet reached the post-core-collapse phase (Trager et al. 


1995). This is to be expected for dynamically young clusters 


with their current half-mass relaxation times, £ r h, of 2.63 
Gyr and 10.47 Gyr respectively ( Harrisp OlO). Our simula¬ 
tions also show that our best-fit models had not undergone 
core-collapse yet. And also for G1 we find that the best-fit 
conditions do not undergo core-collapse within its evolution 
up to its observed age. As argued in Section [3j given Gl’s 
half-mass relaxation time, it is reasonable to assume that 


G1 has not undergone core-collapse yet (Baumgardt et al 


2003b). 


4.3.4 Accuracies on observables 


In Figure 18 we study how including smaller and/or larger 
errors on the observed mass and half-mass radius will effect 
the determined initial conditions. The figure shows the 
initial and final condition distributions in half-mass radius 
versus total mass for the cluster M22 for three different 
error percentages. We chose M22 for its characteristic 
initial condition morphology consisting of shapes (i) and 
(iii) and we wanted to see if this morphology is conserved 
if the observables had larger or smaller errors. Decreas¬ 
ing/increasing the errors also decreases/increases the size 
of the 99.7% confidence region, but the degeneracy in the 
initial half-mass radius and the overall morphology remain 
the same. This makes it more likely that the result that 
M22 could both be a core-collapsed cluster and not core- 
collapsed cluster is robust and that it is not a consequence 
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Figure 18. This figure is similar to the lower panels of Figure[9] but this time we show the 99.7% confidence contours for the simulation 
fitted to observations without a correction for mass segregation, using a 5% (black), 10% (red) and 20% (blue) error on the observables. 


of not being able to distinguish between the two due to 
observational errors. We furthermore see that the lower the 
observational errors in the parameters mass and half-mass 
radius, the better one can especially constrain the initial 
mass of that cluster, but also slighty the half-mass radius. 

In Figure [19] we compare our 2D simulations with our 
5D simulations to infer the stability of the probable initial 
conditions against observational errors in the galactocentric 
radius, orbital velocity and age - in essence combining the 
results from figures like Figures [6] and [15] per cluster. We 
can see that including distance, age and orbital velocity as 
nuisance parameters broadens the 99.7% confidence regions 
for each of the clusters, and that this broadening is most 
visible for the clusters with a degeneracy (morphologies (i) 
or (ii)) in the high probability part of their initial conditions. 
This is because these clusters undergo either expansion or 
contraction and especially these processes are sensitive to 
the parameters setting the tidal field (Rgc and v) and the 
time the time that it could have expanded/contracted (r). 
The morphology is conserved for all but one cluster; in 
these cases the 99.7% regions are just broadened upward 
and downward. Upward if the variation in one of the 
nuisance parameters increases the strength of the tidal field 
(smaller Rgc and larger v); in these cases the clusters must 
have started out out more massive to withstand the larger 
amount of mass loss. An upward shift can also be caused 
if the sampled age is larger than r Q bs- This is because the 
cluster spends a longer time losing mass and thus must 
have started out with a higher mass to reproduce the same 
observed mass. In the opposite cases (decreasing the tidal 
field strength or smaller ages) we see a downward shift. See 
especially the clusters M67, NGC 6397, M4, Pal 14 and 47 
Tuc in Figure [19] M22 also has a degeneracy towards small 


initial half-mass radii, but the effects mentioned above are 
much weaker here. This could be due to the lower sampling 
of initial radii log(rhm/pc) < 0.0. 

We also observe a less prominent, but still visible shift 
in the initial conditions of the degeneracy-free clusters 
Pal 4 and G1 in Figure [19] There we see that a stronger 
tidal field, respectively a shorter evolution time, causes 
the initial condtions distribution to shift to slightly larger 
half-mass radii, indicating that slightly larger clusters are 
more probable to reproduce the observables in these cases. 
For the cluster uj Cen this shift due to the increase of 
tidal held strength even goes that far, that the morphology 
changes: (iii) —> (ii) & (iii). This indicates that initially 
larger (and more massive) clusters which contract during 
their lifetime can also reproduce the observables for this 
cluster. Ah taken together this shows that for some clusters 
accurate observations (error bars < 10% of the mean of 
the observable) are required for determining the initial 
conditions, whereas for the other clusters smaller accuracies 
will be sufficient. 


4-3.5 Star cluster evolutionary tracks 


For our best-fit models with or without a correction for 


mass segregation (see Section 2.3) we make cluster evolution 
tracks by plotting the mass and half-mass radius at different 
equally spaced time intervals in a mass versus half-mass 
radius diagram, see Figure [20] Just before submitting we 
noticed the recently submitted paper by Pfalz ner et al.| 
(2014) in which the authors also make evolutionary tracks 


for young massive clusters covering the first 10 Myr. These 
authors simulate the formation and subsequent expansion 
due to redisual gas expulsion. The evolution of a cluster in 
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Figure 19. A comparison of the 2D and the 5D results to infer the stability of the initial conditions against observational errors in 
the Galactocentric radius, orbital velocity and age. The black (grey) lines show the 99.7% confidence contours of the initial (final) 
conditions of the 2D results. The red (yellow) lines show the 99.7% confidence contours of the initial (final) conditions when including 
the Galactocentric radius, orbital velocity and age as nuisance parameters, i.e. the 5D results. 


the mass-radius diagram during the first 20 Myr is shown 
in their Figure 4. Considering the fact that we model star 
cluster evolution after all redisual gas has been removed and 
the cluster is back in virial equilibrium, our evolutionary 
tracks would, in principle, start some time after their 
tracks end, assuming it would still take some time for the 
cluster to re-virialise. Another paper where the evolution 
of star clusters is studied by means of evolutionary tracks 
is Kiipper et al. (2008), where the authors construct a 


dynamical temperature-luminosity diagrams and show that 
most of their investigated cluster families share a common 
sequence in this diagram. 


We additionally plot the lines of constant half-mass 
relaxation time and constant half-mass density. We calcu¬ 
late the half-mass relaxation time according to Eq. (17) 
of |Portegies Zwart et al. (2010), using r v — rhm for virial 
radius, r v , and In (A) H 10 for Coulomb logaritm, In (A) 
= ln( 7 iV). This is a reasonable assumption for all of the 
clusters using 7 = 0.02 ( |Giersz &; He ggie 1996) and having 
N vary in the range 10 3 — 10 Y . 


From Figure [20] we see that most best-fit cluster models 
start their evolution with a rapid mass loss and expansion 


phase, due to stellar evolution. After that, most of the 
clusters continue to expand their half-mass radius and lose 
mass, but at a much slower rate. As a logical consequence 
the half-mass density of these clusters decreases, whereas 
their half-mass relaxation time increases. Some clusters 
continue at this pace for the rest of their evolution (Gl, uj 
Cen, 47 Tuc, M22 and Pal 4 in both cases, and Pall4 when 
we did not correct for mass segregation). 


When we did correct for mass segregation, the best-fit 
models for M67 and Pal 14 are solely contracting during 
their life time. The interesting feature here is that both 
clusters seem to contract approximately along a line of con¬ 
stant log(phm)- The clusters exhibit mass loss due to stellar 
evaporation over the tidal radius leading to the decrease 
of the half-mass radius (Henon 1961 Gieles et al. 2011) 
and, as these authors mention, this contraction happens at 
(roughly) constant half-mass density. The best-fit models 
with a correction for mass segregation for NGC 6397 and 
M4 start out more than three orders of magnitude larger, 
but at lower mass, compared to the models without a 
correction for mass segregation. Their evolution also starts 
out with a phase of rapid mass loss due to stellar evolution 
and the associated expansion of the half-mass radius. This 
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Figure 20. The star cluster evolution tracks in a total mass vs half-mass radius diagram for our best-fit model for each of the nine 
validation clusters in the case that we did not correct the observations for mass segregation (left panel) and the case where we did right 
panel. For each cluster the square with error bars shows the observed mass and half-mass radius with 10% errors in both parameters, the 
large solid circle shows the initial condition for our best-fit model and the solid triangle shows the final condition. The solid line shows 
the evolution of the best-fit model in mass and half-mass radius and the small solid circles over-plotted on this line mark a passed Gyr of 
evolution, such that M4 has twelve small over-plotted solid circles and M67 has four. The dotted lines show lines of constant half-mass 
relaxation time and the dashed lines show lines of constant half-mass density, corresponding to the values indicated. The shaded regions, 
the lines depicted by green triangles, cyan stars and turquoise hexagons are as in Figure [5] 


is followed by a tidally-limited contraction phase, with the 
decrease of the half-mass radius while log(phm) remains 
close to 2. 


We see that the best-fit models without a correction 
for mass segregation for M67, NGC 6397 and M4 start out 
being very compact (see the beginning of Section 4.3 for a 
discussion on their high initial densities). This is followed 
by a relatively slower expansion of the half-mass radius 
and a slower mass loss phase. After about 1 Gyr, for M67, 
and 3 Gyr, for NGC 6397 and M4, these clusters enter a 
phase where the amount of expansion decreases and comes 
to a halt. From this point on the track bends towards 
the contraction of the half-mass radius. It seems that the 
models of uj Cen, 47 Tuc and Pal 14 without a correction 
for mass segregation are just about to enter that phase and 
the models for uj Cen, 47 Tuc and M22 with a correction 
for mass segregation. 


4-3.6 Determining initial conditions for observed star 
cluster 

With our emacss-mcmc method we are able to reproduce 
the cluster observables well, i.e. reaching a maximum 
posterior probability p > 0.99, for all of the 5D runs and 
most of the 2D runs. However, the 2D runs for NGC 6397 
and M4 without a correction for mass segregation and 
for Pal 14 with a correction for mass segregation reach 
a maximum posterior probability of ~ 0.34, ~ 0.78 and 
~ 0.46, respectively; see also Figures [7| [8] and [T2| where we 
notice that the observation data points are not included 
in the 99.7% a confidence regions in the the mentioned 
cases. Here we see that for M67 the best-fit cluster model 
evolves to a slightly smaller cluster than observed, and 
for NGC 6397 and M4 to slightly larger clusters, but 
with similar masses as the three observed clusters. For 
Pal 14 our best-fit cluster without a correction for mass 
segregation evolves to a final radius in the observed range, 
but significantly more massive than the observations. But, 
we see similar offsets in final mass and/or final half-mass 
radius for the dedicated study for these clusters, see Table [ 3 ] 


The reason for not being able to reproduce the clus- 
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ter parameters can be indicative of two things: some of 
the modeling assumptions are incorrect or the observables 
are poorly constrained. The first point is usually what is 
going on for modeling star clusters, since the model at 
hand always uses a number of assumptions. The second 
point is also interesting, nevertheless. For example, if the 
assumptions made are proven to be reasonable based on 
the observations, e.g. that the cluster’s orbit is observed 
to be very close to circular, and we know that one of the 
observables, e.g. the orbital velocity, v, is determined with 
a large error bar. If we then run the method and are not 
able to reproduce the cluster observables, one can run the 
method with v as a nuisance parameter, probing different 
values of the orbital velocities within the error bars. When 
the second run is able to reproduce the cluster observables, 
this shows that the orbital velocities for which this better 
match to the observations is found, are more suitable values 
for the observed orbital velocity within the observed errors. 

We see something similar happening for the above 
mentioned clusters, where we are not able to reproduce the 
cluster observables in the two-dimensional run, but are able 
to do so in the five-dimensional runs. However, we can not 
exclude the model assumptions to be responsible to the 
mismatch in the 2D runs. Moreover, we also see that this 
mismatch only occurs for one of the cases only, e.g. M4’s 
observables are not reproduced without a correction for 
mass segregation, but are if we do include the correction. 


and 68.3% confidence regions and thus agree well with 
the dedicated studies. For M4 and NGC 6397 our method 
shows, though, that there is a wide range of possible 
initial radii and some spread in initial mass, which are 
also probable. For 47 Tuc our modeling shows that the 
most probable initial conditions well constrained in a 
two-dimensional gaussian-shaped region, consistent with 
the best-fit model of Giersz & Heggie (201 1|) . For M22 the 


best-fit initial condition of Heggie & Giersz ([2014 ) is on the 


2011). Foi 
ersz ( |2014 


high initial half-mass radius border of our 99.7% confidence 
region, but in the same mass range. Our methods thus 
favors smaller initial clusters for M22. For cu Cen the best-fit 


initial condition of Giersz & Heggie (2003) is not included 
in our 99.7% confidence regions. Even though the mass 
range is comparable, our method shows a favor for smaller 
initial half-mass radii. The final mass of |Giersz~fc Heggie 
(2003) and our final mass are in the same range. 


For Pal 4 three out of seven initial conditions of the 


model we compare to of Zonoozi et al. (2014) are included 
in our 99.7% confidence region. This shows that both 
models agree in the initial mass of the cluster, but our 
modeling favors smaller initial half-mass radii. However, 
the final conditions of all seven initial conditions are within 
the 99.7% confidence region of the final conditions. As 
seen in Section [4.2[ this seems to indicate that the N-body 
modeling of Pal 4 led to less expansion than the EMACSS 
modeling. 


4-3.7 Comparison to dedicated studies 

We now determine how probable the best-fit initial con¬ 
ditions of the dedicated studies are according to our 
independent emacss-mcmc results, by checking whether 
these initial conditions are included in our 99.7% and/or 
68.3% confidence regions for the simulations fitted to 
observations without a correction for mass segregation, see 


For Pal 14 a fraction of the initial conditions of [Zc 


3 


Section 2.3) 


For M67 the best-fit initial condition of Hurl ey et al.| 
(2005) is clearly outside our confidence regions, but the 


differences and similarities of the models are not univocal. 
Our modeling shows a degeneracy, more significant in the 
initial half-mass radius, but also in the initial total mass. 
Therefore, when we consider our models with the same 
initial half-mass radius as the best-fit N-body model, our 
models are significantly less massive. When we consider 
the similar initial total mass, our models are either more 
compact or more extended. However, this wide range of 
initial conditions allowed by our modeling all evolve to a 
small region in half-mass radius, but with a spread in final 
mass. The final condition of Hurley et al. (2005) is found 


in this confidence region as well. Compared to our least 
massive models, this again shows the more prominent mass 
loss in the iV-body simulation compared to the evolution 
with emacss. But compared to our more extended initial 
clusters, it shows a larger amount of expansion, since it 
started out with a smaller initial half-mass radius, but end 
up with similar final half-mass radii. 

For NGC 6397, M4 and 47 Tuc the best-fit initial 
and final conditions of GH03-14 are included in our 99.7% 


et al. (2011) are included in the high mass tail of our 99.7% 


confidence regions. In this case both models agree on the 
initial half-mass radius, which shows to be degenerate in 
our modeling, but our method favors a lower initial mass 
and allows even more extended initial conditions. Since 


most of the final conditions of Zonoozi et al. (2011) agree 


quite well with our final conditions, albeit on the high-mass 
end, this again shows what seems to be a general trend: the 
clusters modeled by direct N-body integration lose more 
mass and expand less. 

For G1 the best-fit initial conditions of Baumgardt 


et al. 


(2003b) are clearly outside our 99.7% confidence 
regions; both their initial mass and half-mass radius is 
smaller than those in of our confidence regions with or 
without mass segregation. This is consistent with the results 
from Section \4 .21 Our model favors more massive and more 
extended initial conditions than the best-fit (scaled) initial 
condition of Baumgardt et al. (2003b). 


4.4 Validation 

One could wonder how useful it is to constrain the initial 


conditions based on only two parameters (Heggie & Giersz 


2008). The reason is that constraining the cluster initial 


conditions based on just mass and half-mass radius is in 
principle not enough to truly pin-point the initial conditions 
of an observed cluster, which we also show with the de¬ 


generate shapes (i)a and (i)b in Section 4.3.3 Additionally 
good fits to the surface brightness, the velocity dispersion 
and the luminosity profiles (at a few different radii) are 
required (Heggie & Giersz 2008]). But if one wishes to 
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obtain a reasonable set of probable initial conditions that, 
given the assumptions, contain the initial condition of a 
particular observed cluster with 99.7% confidence, to have 
a good starting point for follow-up modeling and if one 
wants to get a first understanding this cluster must have 
evolved and how that depends on several input parameters, 
our method provides a decent way to do this. 

We see that the emacss-mcmc method does a remarkably 
satisfying job in finding initial conditions for observed star 
clusters. In the direct comparison we found good agreement 
for most of the clusters. For the two clusters uj Cen and M67 
we found poor respectively no agreement, but we argue that 
this comes from the differences in the underlying physics 
of emacss and the codes used in these dedicated studies. 
In the independent emacss-mcmc runs we were able to 
evaluate whether the best-fit initial conditions found by the 
dedicated studies were also good according to our method 
or whether they could be improved. And that is the main 
strength of our method: being able to evolve a distribution 
of initial conditions, study degeneracies and get a good 
grasp on which set of initial conditions are appropriate for 
a given observed cluster, such that these can be followed 
up with more detailed modeling, such as Monte Carlo or 
TV-body simulations. 

Moreover, we have shown with our method that we 
are able to add more dimensions to the initial conditions 
space while preserving the performance and speed of the 
simulation. We have shown this for a five-dimensional 
initial condition space to sample from, see Section [4.3| but 
one could obviously add more parameters to the initial 
condition’s space to fit for, but one could also add nuisance 
parameters over which one marginalizes; this something we 
want to explore more in future work. 


5 SUMMARY 

In this paper we presented our emacss-mcmc method 
with which we are able to derive a distribution of initial 
conditions that, after evolution up to the cluster’s current 
age, evolves to the currently observed conditions. We 
validate our method by applying it to a set of star clusters 
which have been studied in detail numerically with iV-body 
simulations and Monte Carlo methods (hereafter: the 
dedicated studies): the Galactic globular clusters M4, 47 
Tuc, NGC 6397, M22, uj Centauri, Palomar 14 and Palomar 
4, the Galactic open cluster M67, and the M31 globular 
cluster Gl. As in the dedicated studies, our derived initial 
conditions are the conditions of a star cluster after residual 
gas expulsion and re-virialisation. We conclude that the 
results of our method are in agreement with the dedicated 
studies for the majority of clusters. For the two clusters 
uj Cen and M67 we find little agreement, but we argue 
that this due to the differences in the underlying physics of 
emacss and the codes used in these dedicated studies. For 
example not having a direct treatment for binaries in the 
parametrized code, which is especially important for the 
evolution of an open cluster. 

We have furthermore discussed the following points: 


• We have shown the importance of correcting the ob¬ 
servations for mass segregation. If one does not correct the 
observed radius for mass segregation and derives initial con¬ 
ditions that, after evolution up to the clusters age, match the 
uncorrected, and thus smaller, half-mass radius, this mod¬ 
eled cluster does not necessarily resemble the actual cluster 
of interest. 

• We have shown that the distribution of initial condi¬ 
tions can contain two types of degeneracies: a) a degener¬ 
acy towards smaller half-mass radii associated with core¬ 
collapse, and b) a degeneracy towards larger half-mass radii, 
associated with contraction. 

• We made star cluster evolutionary tracks for our best-fit 
models and discussed how the different phases of the cluster 
evolution are distinguishable in these tracks. 

• We find that there is a connection between the mor¬ 
phology of 99.7% confidence region of initial conditions and 
the dynamical age of a cluster and that a degeneracy in 
the initial half-mass radius towards small radii is present for 
clusters which have undergone a core-collapse during their 
evolution time. 

We conclude that our emacss-mcmc method does a satisfy¬ 
ing job in finding initial conditions for observed star clusters 
which can be used in follow up modeling. In forthcoming 
work are applying our method to two groups of star clusters 
in order to provide initial conditions which can be followed- 
up with more accurate methods: Galactic and extragalactic 
globular clusters. 
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Figure Al. The same as Figure p~5] but for the cluster NGC6397. 
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Figure A2. The same as Figure [l5| but for the cluster M4. 
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Figure A3. The same as Figure [6] but for the cluster M22. 
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Figure A4. The same as Figure [T5| but for the cluster Pal 4. 
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Figure A5. The same as Figure [T5| but for the cluster 47 Tuc. 
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Figure A6. The same as Figure pL5| but for the cluster Pal 14. 
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Figure AT. The same as Figure pL5| but for the cluster uj Cen. 
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Figure A8. The same as Figure [T5] but for the cluster Gl. 














































































